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ABSTRACT 


This report is concerned with incremental methods for computer 
graphics. The application of the incremental approach to some advanced 
problems in computer graphics is discussed and demonstrated. In 
Section II, the problem of fast curve generation and display is discussed. 
This section uses the mathematical approach to curves as a perspective 
projection of polynomial-curves in higher-dimensional spaces. Section 
III, also discusses a fast generation of curves, but using another mathe¬ 
matical approach, the linear-differences method, only two-dimensional 
curves are discussed. Section IV, shows how to make the Warnock 
algorithm, for hidden lines elimination, incremental. Section V, dis¬ 
cusses the generation of half-tone images, in real-time. The technique 
suggested there requires a special-purpose hardware to be built. Section 
VI, discusses an incremental method for finding the intersection of given 
line-segments. The technique suggested there also requires a special 
hardware for implementation. 



TABLE OF CONTENTS 


p afie 


SECTION I : INTRODUCTION I 

SECTION II : DISPLAY OF CURVES 2 

II. 1 Introduction 2 

II. 2 On 2D Curves (Perspective Approach) 3 

II. 3 On 3D Curves (Perspective Approach) 5 

II. 4 The Incremental Method for Fast Generation of 

Polynomials 5 

II. 5 The Errors in the Incremental Computation of 

Polynomials 12 

II. 6 Comparison of the Perspective and the Linear 

Differences Methods for Curve Generation 15 

II. 7 Mathematical Justifications 18 

SECTION in : LINEAR DIFFERENCES CURVES 25 

III. 1 Summary 25 

III. 2 The Equivalence of the Various Definitions 27 

III. 3 Obtaining the Generating Matrix for a Conic (Implicit 

Form) 29 

III. 4 Obtaining the Generating Matrix for a Conic 

(Geometric Approach) 31 

III. 5 Characterization of Conics by the Generating Matrices 33 

III. 6 On the Spacing Between Successive Points 36 

III. 7 Straight Lines as Degenerate Conics 38 

III. 8 On Curves Created by Matrices with a Non-unit 

Determinant 39 

III. 9 Programming Aspects of Coding the Linear Difference 

Scheme 41 

III. 10 Examples 42 

SECTION IV : INCREMENTAL METHODS FOR HIDDEN-LINE 

ELIMINATION 56 

IV. 1 Introduction 56 

IV. 2 A Brief Discussion of the Warnock Algorithm 57 

IV. 3 The Incremental Approach to the Warnock Algorithm 58 

IV. 4 The Incremental Solution for the In/Out Problem for 

Polygons 63 


v 




Page 

SECTION V : PRODUCTION OF HALF-TONE IMAGES IN 

REAL TIME 70 

V. 1 Introduction 70 

V. 2 Line Controller for Generalized Coordinates 71 

V. 3 Production of TV-scan Images 73 

V. 4 Production of Radar Scan Images 78 

SECTION VI : FAST INTERSECTING OF LINES 83 

VI. 1 Introduction 83 

VI. 2 The Sutherland Interpolator 83 

VI. 3 Lines Intersector 90 

VI. 4 On String-Line Intersecting 97 

REFERENCES 103 


vi 



LIST OF FIGURES 


Figure No. 


Pa ge 

II. 4. 1 

Schematic drawing of the hardware for the 
forward differences technique 

9 

II. 4. 2 

The behavior graph of the system in Fig. II. 4. 1 

10 

II. 4. 3 

The hardware for the backward differences 
technique 

11 

II. 4. 4 

The behavior graph of the system in Fig. II. 4. 3 

11 

III. 4. 1 

2 2 

The ellipse ^ = 1 

\ M. 

32 

III. 4. 2 

A tilted hyperbola 

33 

III. 10. 1 

A circle (generated by the LDM) 

44 

III. 10. 2 

A circle (generated by the PM) 

45 

III. 10. 3 

Two ellipses 

46 

III. 10. 4 

A family of hyperbolas 

47 

III. 10. 5 

A circular spiral 

48 

III. 10. 6 

An elliptic spiral 

49 

III. 10. 7 

A star 

50 

III. 10. 8 

A star spiral 

51 

III. 10. 9 

A quadratic parabola 

52 

III. 10. 10 

A cubic 

53 

in. io. n 

An oscillated hyperbola 

54 

III. 10. 12 

A family of curves 

55 

IV. 3. 1 

A window subdivided into 4 subwindows 

58 

IV. 3. 2 

Subwindows divided into subwindows 

59 

IV. 3. 3 

A tree representation of Figure IV. 3. 2 

60 

IV. 3. 4 

The relations between the polygon/window relations 62 


Vll 





Figure No. Page 

IV. 3. 5 The simplest Warnock algorithm for line drawing 64 

IV. 4. 1 A polygon describing the digit n 4 n 66 

IV. 4. 2 A Karnaugh map for the signs conditions 67 

IV. 4. 3 Examples for window/line relations 69 

V. 2. 1 Intensity profile 72 

V. 2. 2 A shader for constants intensities 74 

V. 2. 3 A shader for linear intensities 75 

V. 4. 1 A scheme for generating X and Y 79 

V. 4. 2 A multiplication free scheme for computing X and Y 80 

V. 4. 3 An {r , </>) shading system 82 

VI. 2.1 A line intersecting the y-axis 84 

VI. 2.2 X scaling does not change the intersection 86 

VI. 2. 3 Using the Sutherland interpolator as a divider 87 

VI. 2. 4 A schematic drawing of the Sutherland interpolator 88 

VI. 2. 5 The states of the Sutherland interpolator 90 

VI. 2. 6 The control logic for the Sutherland interpolator 91 

VI. 3. 1 Intersection of segments, which are not parallel 

to the axes 93 

VI. 3. 2 Figure VI. 3. 1 after the rotation 93 

VI. 3. 3 The segments {P^, P^} are trivially rejected 96 

VI. 4. 1 A line and a string 98 

VI. 4. 2 The line/segments intersector 100 

VI. 4. 3 The operation of the line/segments intersector 101 

viii 




SECTION I 


INTRODUCTION 

The spirit of the incremental methods is to organize the steps 
of a computation so that each step can use information prepared by 
earlier steps. Such procedures eliminate redundant repetitive 
computation, and simplify each computation step. 

In this thesis I will exhibit five examples of incremental methods. 

The examples shown here all apply to computer graphics. Computer 
graphics is particularly responsive to the application of incremental 
methods because a large amount of information must be processed 
quickly to produce a picture, while the basic operations to be per¬ 
formed on the data are relatively simple. Computer graphics is 
also particularly receptive to incremental techniques because con¬ 
ventional general-purpose computers do not do this job nearly as 
well as it can be done by incremental techniques. Moreover, the 
incremental approach is applicable to a wide variety of computing 
problems, and not to computer graphics only. 

The five examples exhibited here are: linear-difference scheme 
for curve generation, perspective scheme for curve generation, the 
Warnock algorithm for hidden-line elimination, half-tone image pro¬ 
duction and a fast lines-intersecting technique. 

The incremental method for perspective curves, which is described 
in Section II, was implemented in the three-dimensional-graphic- 
hardware project at Harvard, and was simulated on the PDP-1 computer. 
The incremental method for generating curves, which is described in 
Section III, is used by a PDP-1 system for drawing curves on a CRT, 
from which the pictures appended.to that section were taken. The 
search-ordering strategy which makes the University of Utah hidden- 
line elimination algorithm incremental was incorporated into their 
hidden-line elimination system by its author, John Warnock. This 
approach is described in Section IV. A'variation of the incremental 
method for half-tone image production, which is described in Section V, 
was implemented in hardware by the General Electric Company. 
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SECTION II 


DISPLAY OF CURVES 


(II. 1) Introduction 

The ability to display straight lines on a CRT is a standard 
feature of these devices. However, many applications require 
curvilinear display. Typical schematic drawings like block-diagrams 
and logic-designs can be drawn with straight lines only, but attempts 
to display real world objects like machined components require curve 
display since the world is not rectilinear. 

Any curved segment can be displayed accurately up to the scope 
resolution by displaying a set of points close enough together on the 
curve. Curve display by separate calculation and separate storage 
of each point is not feasible for many applications which are restricted 
by real time requirements or by storage limitation. A similar diffi¬ 
culty, existing in the past for straight lines, has been resolved by the 
introduction of line-generators which generate all the points of a line 
segment^ given its end points. Carrying this approach forward to 
curves leads to a "curve-generator n , based upon a small number of 
parameters, which can generate a curve segment in real time. 

Two incremental curve generators are suggested in this work. 

One is based on a three-dimensional perspective approach (Section II) 
and the other on a two-dimensional approach (Section III). Both 
generate points on the curve so that these points may then be connected 
by a line generator. 

Each of these curve generators can generate a family of curves. 
The two families are not the same, but both include all the conics. A 
comparison of these generators is made in (II. 6). The mathematical 
justification of our methods, and of the assumptions used in this section 
can be found in (II. 7). 


up to the scope resolution 
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(II. 2) On 2D Curves (Perspective Approach) 


As the displaying screen is a 2D surface, 2D curves are the core 
of any curve representation. The ability to specify and display curves 
in 2D is the key to representing curves in any space. In (II. 7. 1), we 
prove that any conic segment in 2D is a perspective image of the 
canonic-parabola in 3 d\ or of its images under linear transforma¬ 
tions. Therefore the problem of specifying a conic segment in 2D 
is equivalent to the problem of finding the linear transformation of 
the canonic-parabola which transforms it into a 3D curve whose 3D 
perspective projection is the desired conic. As any linear trans¬ 
formation can be realized by a matrix multiplication, the problem of 
specifying a conic segment is replaced by the problem of finding the 
matrix of the linear transformation. This concept of perspective 
images was already discussed in the late 19th century^’ 


The family of conics includes ellipses and circles, hyperbolas, 
parabolas and straight lines as a degenerate case. Conics never 
intersect themselves and never have an inflection point, as shown in 
(II. 7. 1). If more general curves are required, one can use an 
"mth-degree-conic 11 which is defined to be the perspective projection 
of a 3D curve whose components are polynomials of degree m in some 
parameter t. "Conics" usually mean "2nd-degree-conics " but I will 
not distinguish them from the general case, in spite of the importance 
of second degree conics for many applications. Properties of the 
mth-degree-conics can be found in an unpublished memorandum by 
Professor S. A. Coons and in [4]. In general, the higher m is, the 
more general is the mth-degree-conic, but of course, more conditions 
are required to specify it uniquely. In order to simplify the notation 
hereafter, we will use m = 3. The changes required for other values 
of m are self-evident. 


The canonic parabola in 3D is the curve given by (x, y, z) = (t , t, 1). 
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The explicit form of the mth-degree-conic (for m = 3) is: 



which can also be written as: 


r 3 y 3 


z 


3 \ 


(x, y, z) = (t 3 ,t 2 ,t, 1) 



z 


1 


z 


0 


/ 


A computational implementation of this representation, for generating 



and using a matrix multiplier for getting the (x, y, z) values. This is 
very useful for devices which happen to have a matrix multiplier 
capable of performing the multiplication of a [lx(m + l)] vector by a 


[(m+l) x3] matrix. The scheme can be improved by generating the 

3 2 

sequence {t , t , t, l} instead of storing it. This sequence can be 


generated very fast incrementally. However generating the polyno 



sequence. The curve generator should therefore generate successive 

values of x(t), y(t), z(t) and project these points on the scope plane by 

2 

dividing x(t) and y(t) by z(t) to get the scope coordinates of the points . 


An £xk matrix multiplier is a device which can store an fxk matrix A 
and multiply it by any given lxf vector V, to find the lxk vector VA. 
This operation requires fxk multiplications and (f-l)xk additions. By 
using parallel processing, a matrix multiplier can execute the multi¬ 
plication in t multiply-times only. Such a 4x4 matrix multiplier is de¬ 
scribed in the final report of Harvard contract XG-2972. 

2 

It is easy to see that if the viewing plane is z=l, and the projection 
center is the origin, then the point (x, y, z) is projected on the point 
x v 

(—, x , 1) on the viewing plane, 
z z 
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The scope coordinates are given to a line generator which generates 
line segments between the curve points while the new points are 
generated. The incremental method for generating successive values 
of {x(t), y(t), z(t)} is described in (II. 4) below. 
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(II. 3) On 3D Curves (Perspective Approach) 

For many applications it is very important to specify curves 
which satisfy conditions in the 3D space. After meeting these condi¬ 
tions, the curves are projected into the 2D space for display. One 
can easily generalize the 2D methods to 3D by introducing one more 
component and treating a 3D curve as a perspective projection of some 
4D curve into 3D space. As in the 2D case, specifying a curve is 
equivalent to finding some matrix. Because of the additional compo¬ 
nent, there are more unknowns than in 2D in the equations defining 
the matrix and more conditions are required to specify curves uniquely. 
The perspective transformation from the 4D space to the 3D space (for 
curve definition) and the transformation from the 3D space to the 2D 
space (for display) can be combined in one perspective transformation, 
in order to simplify the displaying process. 


Consider the following example: let {x(t), y(t), z(t), w(t)} be a 
curve in the 4D space. Its perspective projection into the 3D space is 
the following curve: , 'w(t )^ w ^ ose ^D perspective projection 

is the following curve: 


( x(t) . Jdil. = /xiil y(*) \ 

\w(t) * w(t) * w(t) * w(t )J ^z(t)’z(t)J 

Note that there is no need to evaluate w(t) at all. There is only need 
to evaluate {x(t), y(t), z(t)} and to display {x (t)/z(t), y(t)/z(t)} , exactly 
as in the 2D case, since the 3D curve is displayed in 2D regardless of 
the dimension of the space in which it is defined. 


(II. 4) The Incremental Method for Fast Generation of Polynomials 

For each of the 3 components (x(t), y(t), z(t)} we have to evaluate 
a polynomial of degree m, at some set of values of the parameter 
{t q, tj , 1 2 . . . t^} . For simplicity we choose the range of t to be [0, 1] 
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and T = — , i. e. , equally spaced values of t. The usual techniques 
of using multiplications to evaluate polynomials take too much time. 
The fast way to evaluate these polynomials is the finite differences 
scheme [4] which is an incremental method, that prepares information 
from each point to be used in the generation of the next point. This 
incremental method requires only m addition for evaluating a poly¬ 
nomial of degree m. These additions can be executed simultaneously, 
and thus need only one addition time. The implementation of this 
method as described here requires m adders for each polynomial. If 
there is only one adder for each polynomial, the method can be modi¬ 
fied as shown below. 

We have 3 polynomials x(t), y(t), and z(t) to evaluate. As the 
generating procedure is the same for all of them, we consider here 
one polynomial only, say x(t). 

We define: 6 = , the space between successive values of t. 

In the usual fashion we define the forward difference operator [5]: 

Af(s) = f(s+6) - f(s). If x^ is defined as x^ - x(t-^) = x(i6) then 

A x i = We define further: A^f(s) = | A^ *f(s)] and A^f(s) = f(s). 

Note that as is well known, A m t m - m! <5 m , and this is independent of 

l 

t. From the definition of A f(s) it follows that 


and 


A x. 

l 


A[A J 


£-1 




„<-i 

A x. 


i+1 





J-l A^ 1 4. / 

A x. . t = A x. + A x 
l+l l 


Write the last equation for £ = 1, Z, 3 and get: 


x i+l = x i + 

A x. 

l 

d=i) 

A x. + , = A X. + 

Z 

A x. 

l 

U= 2) 

.2 2 

3 


A X. + 1 = A X. + 

Ax. 

l 

(11=3) 


This can be written as: 


6 


x \ 

A x 

a 


wJ 


i+1 


1 0 o\ / x \ 

0 110 Ax 

0011 A 2 X 

\0 0 0 1 / 'A 3 x/. 


or F , = SF. 
1+1 1 


where F. is the Forward differences vector at x = x., and S is the 

i — i 9 

above matrix. 

We define the backward difference operator by Vf(s) = f(s) - f(s-S). 

Of Of f f 

V and V are defined similarly to A and A . In general V f(s) = Af(s-f<5). 

It is easy to see that ^x. . , = ^x. + V^x. , , . 

3 i+l i l+l 

Substituting f = 1,2, 3 gives: 


X i+1 


x. + V x. . , 
i l+l 


V x* 1 = V x. + V K 
i+l i l+l 

y2 x. , 1 = V 2 x. + V 3 X . , , 
i+l l i+l 


which implies: 


V 2 X . . , = V 2 X . + V 3 X . , , = V 2 X . + V 3 X . 

i+l i i+l l l 


x 


i+l 


x i+l 


= V x. + V^x. . , = V x. + V^x. + V^x. 


i+l 


X. + V X.J ! 

1 1+1 


x. + V x. + V^x. + V^ x 
l l l i 


This can be written as: 



\V 3 x /. 


i+l 


/ x \ 

V x 

V 2 X 

W 3 X / 


or B. = TB. 
1 + 1 1 


where B. is the Backward differences vector at x = x. and T is the 

i — i 

above matrix. We have in mind placing the current values of x., 

m f 1 

Vx. ... V x. (or the V x .) in a set of fast registers. The successive 
i i i 
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values in the register are used for the display and the other regis¬ 
ters are used for conveying information from each point to its 
descendants. The A m x register does not have to be modified as t 
changes over its range, because A m t m = V m t m = m[( 
factorize S and T to indicate the computational procedure: 


i !6 m . Let us 


S = 


T = 


\ 

1 

0 

> 

0 


1 

0 

0 

> 

0 


r 

1 

0 

0 

> 

0 


r 

1 

1 

0 

0 

0 

1 

1 

0 


0 

1 

0 

0 


0 

1 

1 

0 


0 

1 

0 

0 

0 

0 

1 

1 


0 

0 

1 

1 


0 

0 

1 

0 


0 

0 

1 

0 

0 

0 

0 

1 

J 


0 

k 

0 

0 

1 

> 


0 

k 

0 

0 

1 

> 


0 

0 

0 

1 

> 

1 

1 

1 

> 

1 


1 

1 

0 

> 

0 


s 

1 

0 

0 

> 

0 


1 

0 

0 

■v 

0 

0 

1 

1 

1 


0 

1 

0 

0 


0 

1 

1 

0 


0 

1 

0 

0 

0 

0 

1 

1 


0 

0 

1 

0 


0 

0 

1 

0 


0 

0 

1 

1 

0 

0 

0 

1 

J 


0 

k 

0 

0 

1 

J 


0 

0 

0 

1 

J 


0 

0 

0 

1 

J 




AiA^As 


t l 1 

A. is the operation of adding the V x (or the A x) register to the V x 
1 £-1 

(or the A x) register. S is equivalent to executing then A^, 

then A^. T is equivalent to the execution of the same operations, in 

reverse order. There are only three additions involved in either of 

these methods. However, there is a basic difference between them. 

In the forward difference scheme, each f, new M value A^x.^ depends 

only on "old" values A^x., but in the backward difference scheme, each 

"new” value V^x. ,, depends on the current values of the registers, 

"old" V x . and "new" V x. . , . It is therefore possible to execute all 
i i+I 

the additions for the forward difference method scheme in parallel, 
consuming only one addition time for evaluation successive values of 
the polynomials. This however requires a great deal of hardware. 

The backward difference scheme is more economical for sequential 
computation, but consumes m addition-times. Schematic drawings 
for the two methods are shown in Figures II. 4. 1 and II. 4. 3. 

The behavior graphs^ - of these systems are shown in Figures 
II. 4. 2 and II. 4. 4. The behavior graph of the forward difference scheme 


Behavior graphs are sort of "occurrence-graph n [18]. 
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Figure 11.4.1 : Schematic drawing of the hardware for the forward differences 
technique. 
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Figure 11.4.2 ; The behavior graph of the system in figure; II . 4 


1 0 

























V 2 X 

r?3 v 

°3 

i r 

, ^ 


V A 

1 


Figure II.4.3 : The hardware for the backward differences technique. 



Figure II.4.4 : The behavior graph of the system in figure II.4.3 
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has 3 parallel paths which indicate the parallel processing. It is 
possible to introduce some parallelism in the backward difference 
scheme, but this requires more hardware than is needed for the 
forward difference scheme. However the critical operation is the 
division. If the division is slower than 3 additions (as it always is 
in digital systems) there is no sense in providing the extra hardware 
which is required by the forward difference method. In this case, 
both the forward and the backward difference method require on 
divide -time . 

I 

The forward difference scheme is initialized by loading Ax.(O) 
into the registers, and the backward difference scheme is initialized 
by loading V^ x ^(0). 

In (II. 7. 2) we show how the initial differences {V^ x .(0)} and 

i 1 

{ A (0)} can be found. 

(II. 5) The Errors in the Incremental Computation of Polynomials 

The iteration process which is described in the preceding 
subsection may introduce some errors due to the finite precision of 
the machine. We cannot, in general, better the precision, but we 
can modify the iterations to minimize the errors where they are most 
important, at the end-points (e. g. , for curve closure). 

The source of the errors is not the iteration process itself 
(which is multiplication free) but the propagation of the errors in the 
initial values of the differences. 

We can change the initial values of the differences in order to 
make the iterations end as close as possible to the prespecified end¬ 
point. We show now a method which converges to an end-point which 
is within-ty machine resolution points from the given end-point, where 
N is the number of iterations. 

Consider first the forward differences scheme. Define: 
e = [1, 0, 0, 0] and S = I + H, 


12 






1 

0 

°\ 

/ x(n6) \ 


0 

0 

1 

0 


A x(n6) 

H = 





and F(n) = 

? 


0 

0 

0 

1 


A x(n6) , 

\o 

0 

0 

0 / 


\ A^x(n6) / 


We showed that F(n) = S n F(0) and x(n6) = eF(n). The computed end¬ 
point is x(N6) = eS^F(O). Let be the given end-point. Hence, the 
error introduced by the iterations is: 

E = x e - x(N) = x e - [x(0) + NAx(O) + 


A 2 x(0) + A 3 x(0)]. 


Dividing the error by ( ] is a correction for A x(0). The remainder 

V* 3 / /N\ . 2 

of this division divided by ( ^ I is a correction for A x(0), and the 

remainder of the second division divided by N is a correction for Ax(0). 

These 3 corrections reduce the magnitude of the error to be less than 

N. By changing Ax(0) by one resolution unit, the magnitude of the error 

N 

can be reduced not to exceed 


This correction process can be programmed as follows: 

E 


(1) x e - [x(0) + NAx(O) + 


A 2 x(0) + (3 ) A 3 x(0)] 


( 2 ) 


(3) A x(0) + e. 


A^x(O) 


(4) 


a. 


w 


(5) A 2 x(0) + -* A 2 x(0) 


( 6 ) 


lIv M 


N 


-> e. 
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Ax(0) 


(7) Ax(0) + e L 

(8) if [E - (”)e 3 - (”)e 2 - Ne,] 

(9) if [E - ( 3 je 3 - (?) e 2 ‘ Ne l) 

For the backward difference 
works by using: 

x(n6) = eT n B(0) where 


>— then Ax(0) + 1 ) Ax(0) 

< - ^ then Ax(0) - 1 )Ax(0) 

scheme a similar correcting method 

T = I + H + H Z + H 3 


and 

T n = I + nH + [n + ]H 2 + [n + 2 + ("j ]H 3 

If we do not want to change Ax(0) ? and we wish also to get the 
right value for Ax(N) then another correcting scheme can be applied. 
Expand F(N) = S^F(O) and get: 

x(N) = x(0) 4* NAx(O) + 

Ax(N) = Ax(0) + 

A^x(N) = 

A 3 x(N) = 

x(0) and x(N) are given by: 
x(0) = d 

x 

x(N) = a + b + c + d 

X X X X 

Ax(0) and Ax(N) can be found directly by: 

Ax(0) = a 6 3 -f b 6 2 + c 6 

Ax(N) = a (6 3 + 36 2 + 36) + b (6 2 + 26) + c 6 

X X X 


(?) A 2 x(0) ♦ (?) A 3 X(0) 
N A 2 x(0) + (?) A 3 x(0) 
A 2 x(0) + N A 3 x(0) 
A 3 x(0) 


14 


2 3 

Using these values, A x(0) and A x(0) are found from the following 
equations: 

A 2 x(0) + A 3 x(0) = x(N) - x(0) - NAx(O) 

N A 2 x(0) + 

2 3 

Using the values of A x(0) and A x(0) which satisfy (as close as possible) 
these equations reduce the errors in x(N) and Ax(N). For backward 
differences a similar method exists. 

(II. 6) Comparison of the Perspective and the Linear Differences Methods 

for Curve Generation 

The LDM (Linear Difference Method) is described in Section III. 

In this section I have described the PM (Perspective Method) for curve 
generation. 

I compare the two methods, the LDM and the PM (for m = 2) accord¬ 
ing to: 

- complexity of points generation, 

- variety of curves and speed of generation, 

- complexity of the mathematics involved, 

- sensitivity to errors. 

The complexity of points generation 

The PM requires 6 additions and 2 divisions per point. If enough 
hardware is available then all the additions can be carried out in parallel, 
and so can the divisions. Hence the time which is consumed by each 
point is between 1-add-time plus 1-divide-time and 6-add-times plus 
2-divide-time s. 

The LDM requires 4 multiplications and 2 additions. If the hard¬ 
ware is available then the multiplications can be performed in parallel, 
and so can the additions. Hence the time which is consumed per point 
is between 1-add-time plus 1-multiply-time and 2-add-times plus 
4-multiply-times. 


A^x(O) = 


Ax(n) - Ax(0) 






If the curves are generated by software it is preferred to use 
the LDM. If hardware is used then the comparison depends on the 
speed and the cost of the components involved. 

Variety of curves and speed of generation 

The PM can generate only conic segments; ellipses, parabolas, 
hyperbolas and straight lines. It cannot generate complete ellipses 
and circles. The LDM can generate complete ellipses and arcs of 
hyperbolas and generalized parabolas (y = x a ) which do not pass through 
the origin or through infinity. In addition the LDM can also generate 
straight lines, elliptic spirals, stars and various other shapes as 
described and illustrated in Section III. 

The LDM generates the conics always at a "good-speed . The 
PM is not always able to generate a conic section at a "good-speed". 
For example, it is impossible for the PM to generate a circular arc 
at a uniform speed. Because of its good-speed-property the LDM 
needs fewer iterations to display a curve than the PM. (See the figures 
in III. 10). 

Complexity of the mathematics involved 

There is a simple method for finding the parametric matrix for 
any conic segment given by its end-points (see II. 7 and [6]). This 
method can be used for finding the matrices which together represent 
a complete ellipse. However, although it is relatively simple to find 
the parametric matrix, as required by the PM, it is not always easy 
to perform the shape invariant transformations [4] which are required 
to improve the generation speed. 

There is no simple way to find the generating matrix for a given 
conic segment, as required by the LDM, because finding the matrix is 
equivalent to finding sines and cosines. However it is relatively very 
easy to find the generating matrix for a conic (not only a segment) from 
its implicit form or from its geometrical properties. 


A curve is said to be generated at a "good-speed M if the distance 
between successive points decreases when the radius of curvature 
decreases. 
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The simple manipulation of curves, such as rotation, translation, 
stretching and scaling do not have to be performed on each data point. 

In the PM all these operations can be performed on the parametric 
matrix. In the LDM only the stretching and rotation (in some cases) 
have to be performed on the matrix. The other operations do not 
need any operation on the matrix, and are carried out by the specifi¬ 
cation of the initial point (and initial difference), because the same 
generating matrix generates a whole family of curves. 

Sensitivity to errors 

Both, the PM and LDM may introduce roundoff errors. These 

roundoffs propagate during the iterations, and might result in a large 

error after N iterations. We show in this section that the error in 

evaluation each polynomial for the PM can be reduced to be less than 
N 

-y- in magnitude at the end points. We can do that because the iteration 
involve only additions which do not introduce new roundoff errors. The 
LDM uses multiplications which are rounded off, introducing possibly 
new truncation errors, of one machine-resolution unit, at each step. 
These roundoff errors are accumulated and multiplied by T n . Hence 
the LDM is more sensitive to errors, then the PM. In some cases the 
LDM becomes so sensitive that it cannot be used. For example: A 
definition of a hyperbola by its implicit form and a point which is very 
close to an asymptote. The reason for this kind of sensitivity is that 
the same generating matrix generates a family of hyperbolas which all 
converge to the same asymptotes. 

The PM does not have a 1-1 correspondence between curves and 
matrices, unlike the LDM. It is possible to change the representing 
matrix of a curve (even by splitting the curve into some sections) in 
order to make the iteration less sensitive. By this operation the sensi¬ 
tivities due to a small denominator, or to small numerators may be 
improved. 

Unfortunately the generating matrices for the LDM are unique 
(up to the generation speed) and do not have any freedom which can be 
used for improving the sensitivity to errors. 
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(II. 7) Mathematical Justifications 

The purpose of this subsection is to justify mathematically some 
results which are used in the section without proof. 

(II. 7. 1) We will define a family of curves, cr , which con- 

3 ’ m, n’ 

n 7 

sists of curves in R , whose components are polynomials of degree 
m in some parameter t. The perspective projection of cr 


into 


R 


n-1 . 


is called here 77 


We will prove the following results: 

111^ xl 

(a) If S e 77^ ^ then S is a conic segment* 

(b) If S is a conic segment then S e 77^ 

(c) There exists S e cr such that for any P € 77 ? there 
exists a linear transformation T, such that P is the 
perspective projection of TS; 

(d) Conics do not intersect themselves. 

Let cr^ be the family of curve segments in R n , with each component 

being a polynomial in t, with degree not exceeding m. Sea - implies 

m . ’ 

S = {s(t)} = {x-, (t), x 9 (t). . . x (t)} where x.(t) = 2 a. .t^. Let 77 be 
1' ° 2 W n' ' V ' j~o ij m, n 

the perspective projection of cr in P n , the perspective space obtained 
n 

from R by dividing each component by x^. Note that P is isomor¬ 
phic to the closed R n 


T-,n . 

P is: 


If S = (s(t)} = {x^(t), x 2 (t). . . x(t)} e <r^ then its projection in 

Xj(t) x 2 (t) x n-1 (t) 


P = (p(t)} = 


X (t) 9 X (t) * X (t) 

n n n 


e 77 


(a) We want to show that if S e 772 3 t * ien S is a conic segment. 
Proof (following L. G. Roberts [6]): Consider 772 3’ the 

following notation: x = Xp y = X 2 and w = x^. Let 

2 

x = x(t) = X2t + x^t + Xq 
2 

y = y(t) = y 2 t + yjt + y Q 


W 


- w(t) = W2t Iw^t + 


w. 
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In short, this can be written as: 


p = p(t) = (x, y, w) = (t ,t, 1) 


where T = (t , t, 1) and A is the above matrix. 
If A is a nonsingular matrix, define: 


x 2 

yz 

w 

X 1 

n 

w 

x 0 

yo 

w 


2 ' 

1 

0/ 


= T A, 


C = A 


-1 


0 


\l 


0 

■2 

0 


l\ 

0 

°/ 


* -1 -1 * -1 
A a A MA 


then: 


p C p = ( t A)C(tA) V = (tA)(A" 1 MA' : " 1 )(A'V ) = T 


0 

0 


0 l\ 

•2 0 


* 

T =0 


\i 0 0/ 

for all t. This proves that the arc tA is a conic segment. If A is 
singular, then there exists a non-zero vector V such that AV = 0, and 

h 

tAV = pV = (x, y, w) b = ax + by + cw = 0. Hence tA is a straight 
line, which is a degenerate conic. 


(b) We proved above that if P € 7T^ ^ then Pisa conic segment (or a line 

segment, which is a degenerate case of a conic). The converse is more 
important: any conic segment belongs to 77 ^ We prove this as follows: 

❖ 

Proof : Let vCv = 0 be the implicit equation of the conic, and let Vq and 
v ^ be the start and end points of the segment. Let v^ be the intersection 

of the two tangents to the conic through v n and v, . Note that the tangents 

2 u 1 3 

do not necessarily intersect in R , but they must intersect in P . 


Consider the matrix 



1 -2 1 \ 


*1 \ 

A = J V = 

0 2-2 


fv^. 


0 0 1 i 


v„ 


/ 

0 / 
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where f is a scalar. It satisfies the following conditions: 

(i) for t = 1, rA = V 1 : (1, 1, 1) A = (1, 1, 1) J V = (1,0, 0) V = vj 

(ii) for t = 0, rA = V Q : (0, 0, 1) A = (0, 0, 1) J V = (0, 0, 1) V = v Q . 

❖ 

(iii) tA is a conic arc because (tA) C (tA) = 0 for all t. 

This is proved by: 


M 




f v j C v j ^ 

fv^Cv^* 

v i c V ' 

fv T 

C 

fv T 

= 

fv T Cv 1 * 

fv^Cv^,* 

fv T Cv 0* 

l v o i 


l v o 1 


l v o Cv i* 

f v Cv ^ 

0 T 

v o c V j 


0 0 a\ 

0 b 0 


y a 0 0^ 

2 

where a = v^Cv^* and b = f v^Cv^,* . 

v lCvi* = v q^ v q >:c = 0 because and are on the conic. 

v^Cv^* = VjCv^ = 0 because is on the tangent to the conic 

through . 

V T^ V 0* = V 0^ V T* = ® because of a similar reason. 

With no loss of generality we can assume that C is a positive definite 

2 

form. In this case, a = VjCv^ = VgCv^* < 0 and b = f v^Cv^* > 0. 
Next we consider the product A C A* = J V C V' 4 ' J* : 



1 -2 1^ 


1 0 0 a\ 


1 0 0^ 


c “c a ' 

J(V C V*)J* = 

0 2-2 


0 b 0 


-2 2 0 

= 

-c 4b 0 


0 l l 


^ a 0 0 1 


1 - 2 '1 


^ a 0 0 1 


where c = 2a + 4b. Hence 

(tA) C (tA)* = t (J V C V* J*)t* = ct 4 - 2ct 3 + ct 2 = ct 2 (l-t) 2 . 
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This expression vanishes for all values of t if c = 0. This con¬ 
dition can be satisfied by the proper choice of f, namely 


~ 2v t Cv t * 

We showed in (iii) that tA is a part of the conic VCV* = 0, and 
we showed in (i) and (ii) that the end points of tA are v^ and v^. 

This completes the proof of the converse claim. 

(c) Corollary : Any conic segment can be obtained from any non¬ 
degenerate conic segment by a linear transformation. 

Proof : Let v^ = rA^ be a non-degenerate conic segment. Let 
V£ = tA^ be another segment, v^ is obtained from v^ by the linear 
transformation T = as v^T = rA^(A^A^) = tA^ = v^. Note 

that if v^ is the "canonical-parabola” v = (t , t, 1) then A^ = I, the 
unit matrix, and T = A^, which means that V£ is obtained from the 
canonic-parabola byusing the representation matrix of v^ as the 
transformation. 

(II. 7. 2) The forward and the backward difference schemes need 

£ £ 

initialization by loading the values of Ax.(0) or V x.(0) to the appro- 

£ 1 1 £ 

priate registers. We will find first A xj(t) then we will find Ax.(0) by 

substituting t = 0. Later we will find V^x.(t) and substitute t = 0 to 
o 1 

get the V x^(0). 

We use the fact that the V and the A are linear operators in the 
evaluation of the differences of the polynomial x(t): 

A t 3 = 3t 2 6 + 3t6 2 + 6 3 

A t 2 = 2t6 + 6 2 

At =6 

A 2 t 3 = 6t6 2 + 66 3 
A 2 t 2 - 26 2 

A 3 t 3 =66 3 . 
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Then: 


3 2 

Let x(t) = a^t 4- a^t 4- a^t 4- a^. 


A x(t) 

= (a x 6 + a 2 6 2 + a^ 3 ) + t(2a 2 6 + 3a 3 $ 2 ) + 3a 3 6t 2 

A^x(t) 

= (2a 2 § 2 + 6a 3 6 3 ) + 6a 3 6 2 t 

A^x(t) 



Substitute t = 0 and get 


x(0) 

= a o 

A x(0) 

= a l 6 + a 2 62 + a 3 ^ 3 

A Z x(0) 

= 2a 2 S 2 + 6a 3 ^ 3 

A 3 x(0) 

= 6a 3 63 ' 


Find the backward differences by: 


v t 3 = 

3t 2 6 - 3t6 2 + 6 3 

V t 2 = 

2t6 - 6 2 

V t = 

6 

V 2 t 3 = 

6t6 2 - 66 3 

v 2 t 2 = 

26 2 

v 3 t 3 = 

6t 3 . 


3 2 

Again let x(t) = a^t 4- 4- a^t 4- a^. Then: 


V x(t) 

= (a 1 6-a 2 6 2 + a 6 3 ) + (2a 2 6 -3a 3 6 2 ) + 3a 3 6t 2 

V 2 x(t) 

= (2a 2 6 2 - 6a 3 6 3 ) + 6a 3 ^ 2 t 


V 3 x (t) = 6a 3 6 3 
Substitute t = 0 and get: 
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x(0) = a Q 

V x(0) = a 2 6 2 + a 3 

V x(0) = 2a 2 6 2 - 6a 3 § 2 
V^x(O) = 6a 3 ^^ 

To summarize : Let the curve be given by 

/ 


[x(t)j y(t), w(t)] = [t 3 ,t 2 ,t, 1] 


The Forward Differences Scheme : 

We showed in (II. 4) that: 


F(n) = 


x 


x 


X 


X 


w 


w 


w 


w 


\ 


= rA 


/ Xn \ 



i 

0 

°\ r 

A x 

n 


0 

i 

1 

0 



0 

0 

1 

1 

\a 3 x y 

n 


\ 0 

0 

0 

1/ 

Then 

we 

showed that: 


I x(0)\ 
A x(0) 
A 2 x(0) 


= S n F(0) 


n 


/ 


X 


F(0) = 


X 


6 3 + b 6 2 + 


X 


X 




6a 6 3 + 2b 6 2 

X X 


6a 6' 
x 


\ 

1 0 0 0 1\ 

\ 

/ a x\ 


1110 

fi 2 \ 

b 

X 


6 2 0 0 

6 

c 

X 

/ 

\6 0 0 0/ 

\ 1 / 

/ 


Q 


D 


Combine together and get: 

[x(n6), y(n6), w(n6)] = [1,0, 0, 0]S n F(0) = [1,0, 0, 0]S n QDA 
which may be verified by checking that: 

[1,0, 0, 0]S n QD = [(n6) 3 (n6) 2 n 6 1] . 


n 

= QDA 0 

\ol 
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The Backward Differences Scheme: 


We showed that: 


B(n) = 





1 

1 

*\ 

v x 

n 


0 

1 

1 

1 

V 2 X 

n 


0 

0 

1 

1 1 

\V 3 x y 
n 


\o 

0 

0 

1 / 


V x(0) 

V 2 X (0) 

iV^x(O) i 


T n B(0) 


where x = x(n6). Then we showed that 
n 


I 0 0 0 1 \ I 

1-110 
-6 2 0 0 

^ 6 0 0 0 


B(0) = 


J 


\l * 

b 


/ a 




J 




Combine together and get: 

[x(nS) y(n6) w(n 6)] = [1 0 0 0]T n QDA 

which may be verified by checking that: 

[1 0 0 0]T n QD = [(n6) 3 ( n S) 2 (nS) 1] . 


QDA I 0 
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SECTION m 


LINEAR DIFFERENCES CURVES 


III. 1 Summary 

We are interested in the family of curves of the form: 

n = {P(s)|P(s) = T S P(0) J 0 < s < oo} 

where T is a 2x2 real matrix, P(0) is the initial point in 2D space, 
and s is a continuous variable. 

These curves can be displayed by generating the sequence of 
points (P(n)} where n is an integer, and connecting successive points 
by straight lines. The sequence (P(n)} can be generated incrementally 
by using: 

P(n+1) = T P{n) . 

The simplicity of this iteration makes it very attractive for digital 
systems involving either a special hardware or conventional program¬ 
ming. 

There are several different definitions for these linear-differences 
curves. The main ones are: 

(a) P(n+1) = TP(n) or P(n) = T n P(0) . 

(b) AP(n+l) = T^AP(n) where AP{n) = P(n+1) - P(n) . 

(c) AP(n) = T 2 P(n) . 

< d > - T 3 P < n > • 

All these definitions are equivalent. In (III. 2) below we prove this, 
and show the connection between T^ T^, T^ and T. Although we use 
only the first of these definitions we want to point out that, for some 
applications, the other definitions might be more appropriate, (or 
even still other definitions). 

In (III. 3) below we prove that any origin-cente red-conic ^ may be 
generated by the above process. The proof is constructive and gives a 


1 2*2 

An origin-centered-conic is defined by ax + 2bxy + cy = d. 
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"recipe" for getting the generating matrix T, for any conic specified 
by its implicit form. As corollaries from this theorem we get the 
generating matrices for circles and hyperbolas. Each of the generating 
matrices obtained by this method belongs to a one-parameter family 
of matrices, all of which generate the same curves but at different 
speeds. The free parameter can be used to control the generating 
speed, for example, by specifying P(l) on the curve; (with P(0) already 
specified). 

In (III. 4) we show that if two curves are obtained from each 
other by a linear transformation, then their generating matrices are 
similar (in the usual mathematical sense). As a result we have a 
method for constructing the generating matrices for ellipses and 
hyperbolas according to their geometrical properties. 

In (III. 5) we show that the condition det(T) = 1 implies that T 
generates: 

(a) an ellipse if |trace(T)| <2 

(b) a straight line if |trace(T)| = 2 

(c) an hyperbola if |trace(T)| >2 . 

Next we are concerned with the "speed” of the conic generation, where 
the "speed" is defined as the distance between successive computer¬ 
generated points. In (III. 6) we show that unlike the perspective method, 
the LDM makes it possible to generate a circle with a uniform speed 
(i. e. , equally spaced points) and hyperbolas with the right kind of 
speed (i. e. , the speed decreases down as the radius of curvature de¬ 
creases). As a corollary from the equally spaced circle generation 
it follows that it is possible to generate ellipses with the right kind of 
speed. Next we show the connection between the area of successive 
triangles A^ and det(T). In particular we show that the areas of all 
the triangles of a conic are constant. This implies that the spacing 
on a conic is necessarily good. 



is the triangle whose vertices are P(n), 


P(n+1) and the origin. 
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It is important to mention that the family of conics includes 
two kinds of straight lines, the ones which pass through the origin, 
and the ones which do not. In (III. 7) we discuss these lines. 

In (III. 8) we discuss the curves which are generated by matrices 
with a non-unit determinant. We show that for any a, the curve 
y = kx a can be generated by the iteration, as well as elliptic spirals. 

Finally, in (III. 9) we discuss some programming aspects of 
coding the iteration. 

(III. 2) The Equivalence of the Various Definitions 

We will show that the 4 following definitions of linear-differences 
curves, are essentially equivalent, in the sense that they define the 
same family of curves. 


(a) P = TP 

n n i 

(b) AP = T, AP , 

v 7 n 1 n-1 

(c) AP = TP 

n 4 n 


dP 

n 


(d) = T P . 

x dn 3 n 


We will show the equivalence of each definition to (a). 

(III. 2. a) (a) implies (b), i. e. , if a curve belongs to the family 
which is defined by (a), then it also belongs to the family of curves 
which is defined by (b). 




(b) does not imply (a), but: 



. = AP +. . . + AP, 

n i 



= (T n + T n_1 + . • .I)AP q + P Q . 
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Assume that 1 is not an eigenvalue of T, then we can write: 

A p 0 = (T - I)P 0 , 

then 

P n+1 = (Tn + Tn_1 + • • * + I >< T ' DP 0 + P 0 = + (P 0 " ^0> 

which shows that (b) defines the same curves as (a) but they might be 

A 

off-centered by E = - P^. The reason for this possible displacement 

is that (a) requires only an initial P^, but (b) requires initial P^ and 
APq. If AP^ = TPq - P^ = (T - I)Pq then E = 0, and there is no center 
displacement. 

If 1 is an eigenvalue of T, then T generates a straight line, 
which obviously can be generated by (a). 

(III. Z. b) (a) implies (c): 

AP = P ^ - P = (T - I)P = T 9 P . 
n n+1 n n Z n 


(c) implies (a): 


P . ,= P + AP =P +T 9 P =(T 9 +I)P 
n+1 n n n Z n Z n 


(III. 2. c) (a) implies (d): 


P(n) = T n P(0) 

= (l!gT)T n P(0) = (£g T)P(n) = T 3 P(n) . 


(d) implies (a): 

dP(n) _ ~ 
dn 3 


P(n) 


P(n) = e nT 3A ; 

substitute n = 0 and get A = P(0), 

P(n) = (e T 3) n P(0) = T n P(0) . 

Note that if T has non-positive eigenvalues then there does not exist a 
real matrix = Ig T, and (a) does not imply (d). 
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For displaying an off-centered curve, one can generate the 
{P(n)} using (a), and add some displacement to each point, or posi 
tion the initial point P(0), and generate the {AP(n)} using (b). The 
latter scheme saves the addition of the displacement to each point, 
and is, therefore, preferred if a special-purpose hardware is not 
available. 


(III. 3) Obtaining the Generating Matrix for a Conic (Implicit Form) 


Theorem: For any given origin-centered non-degenerate conic, there 
exists a one parameter family of matrices {T(k)} which 
generate the conic, and det[T(k)] = 1 for all k. 


Proof: 


Let the conic be 


P*CP = P* 



P = d . 


We construct a matrix T, such that if P. is on the conic, then so is 

’ 1 ’ 

Define: 


p i+ i 


T P.. Consider P. , , = P. + AP.. 
i l+i i i 


E = (P. + AP.)* C (P. + AP.) - P. C P. 

i i i l i i 


Expand it: 

E = 2ax(Ax) + a(Ax)^ + 2by(Ax) + b(Ay)(Ax) 

2 

+ 2bx( Ay) + b( Ax)( Ay) + 2cy(Ay) + c(Ay) 


= (Ax)- [2ax + 2by + a(Ax) + b(Ay)] + (Ay). [2bx + 2cy + b(Ax) + c(Ay)]. 
For Ax, Ay and any k which satisfy the following: 

Ax = k[2bx + 2cy + b(Ax) + c(Ay)] 

Ay = -k[2ax + 2bx + a(Ax) + b(Ay)] 

E vanishes, which means that P^j is on the conic. Separate Ax and Ay: 

(1 - kb)Ax -kc Ay = 2k[bx + cy] 
ka Ax + (1 + kb)Ay = -2k[ax + by] 
which is in matrix form: 
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Introduce: 


/ 

1 - k 




b c 

-a -b 


(AX\ 

= 2 k 

b 

j 

(x\ 

Uy/ 


i-a 

-bl 

ly 1 


G = 


0 1 

-1 OJ 


a b 1 
b cl 


b c 
-a -b, 


then. 


(I - kG)AP = 2kGP . 


The matrix (I - kG) is invertible for all (except two, at most) values 
of k: 

-1 


AP = 2k(I - kG) GP 


-1 


P. + 1 = P + AP. = [I + 2k(I - kG) G]P. . 

-1 2 

Hence T(k) = I + 2k(I - kG) G. Substitute a, b, c and define e = b - ac. 

Then for / e ^ : 


T(k) = —4- 
1-k e 


l+2kb+k 2 e 


2kc 




l-k 2 e 


1-2kb+k e 

It is easy to verify that det(T) = 1 ? for all k. 
Consider the trace of T(k): 


trace(T) ~ 


trace(T) = 


^1+k^e 

^ 2 
1 -k e 

/ 

< 2 if 
< = 2 if 
> 2 if 


e < 0 
e = 0 ) 
e > 0 


for all k . 


\ / 

It is well-known that e < 0 implies an ellipse, e > 0 implies an hyperbola, 
and e = 0 implies straight lines. 
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2 2 

Corollary : The generating matrix for the circle x + y = a is obtained 
by substituting a = c = 1 and b = 0 (e = -1): 


T = 

l C 


1+k 


/ 1-k 2 2k \ 

I 

\ -2k 1-k 2 j 

\ 


cos 0 -sin 0 

sin 0 cos 0 , 


2 2 

For the hyperbola x - y = (3, substitute a = -c = 1 and b = 0 (e = +1): 


1 

1 1+k 2 -2k \ 

^ ch 

sh 4 ^ 

T = 1 

H i-k 2 

^ -2k 1+k 2 J 

sh jif 

ch<4 j 


-2k -2k 

where 0 = arctg - x and 0 = argth - x . 

1+k 1-k 

(III. 4) Obtaining the Generating Matrix for a Conic (Geometric Approach) 

Theorem : If T generates the curve IT , and the linear transformation H, 
maps [1 into the curve 2, such that 2 = HH then the generating 
matrix of 2 is S = HTH 

Proof : Let II = (P(n)} and 2 = (S(n)} such that S(n) = HP(n) for all n, 


then: 


S(n+1) = HP(n+l) = HTP(n) = HTH _1 S(n) . 


QED 


Namely, S and T are "similar matrices", and have the same eigenvalues 
(and therefore the same determinant and trace). 


This theorem suggests another method for finding the generating 
matrices for ellipses and hyperbola. Consider the ellipse, whose axes 
are parallel to the X-Y axes, the length of the horizontal axis is 2X, 
and the length of the vertical axis is 2p, as shown in Figure III. 4. 1. 
This ellipse is obtained from the unit circle by the transformation: 


l\ 


D = 


\0 



If this ellipse was tilted by the angle a then it could be obtained from 
the unit circle by the transformation R(a)D, where 
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2 2 

Figure III. 4. 1 : The ellipse ^ = 1 


R(ci) = 


cosa -sina 
sina cosa 


-1, 


Hence the generating matrix of the ellipse is E = R(a)DR(0)D R(-a). 

ig matrix of 
^ch^ sh^ 


2 2 

The generating matrix of x - y = 1 is 


t h^> = 


; sh^ ch </> 


Consider the hyperbola in Figure III. 4. 2 ? which is obtained from 
2 2 

x - y =1 by the transformation 


H = R(a) D = 


f cosa 


sina 


- sina 


cosa 



Hence its generating matrix is: 


-1 


T = R(a)DT H (/)D R(-a) 
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Figure III. 4. 2 : A tilted hyperbola 

(III. 5) Characterization of Conics by the Generating Matrices 

Theorem : If n = {P(s) |P(s) = T S P^} and det(T) = 1, then H is a 
conic, whose type depends on trace(T). 

Proof : Assume that T / tl as these cases generate either P^ re- 

peatedly or P^ and -P^ alternately. Let 

’■(! 3 ■ 

Consider the following cases: 

(a) |a+d| > 2 (b) |a+d) = 2 (c) |a+d| < 2 . 

Case (a) : If ja+d ) >2. C onsider X, an eigenvalue of T, 

\ ^ +Y|(a+d) 2 - 1 . 

T can be written as: 


“ 1 1 

i , i ,\ 

f >> r _d ' 

l x 

0 ^ 

T = QDQ where Q = 

U 

oj 

1 

and D = 

\° 

1 

\ / 
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D can be written as 


D = ST^S ^ where S = 


1 V 


and T h = 


ch <f> 

sh <f> 

sh <f> 

ch <f> 


where ch ^ ^ (X + ^ ) and sh ^(\ - Hence, 

T = QDQ -1 = (QS)T h (QS) _1 . 

Tpj generates the hyperbola H, and T generates the hyperbola (QS)H. 
Case (b) : If a+d = 2 and c / 0 then: 

,-1 


a-1 1 \ /1 l\ / a-1 1 


T = 


= QLQ 


0/ 0 1 


The matrix L = 


1 1 


lo 1, 


generates lines since 


P * = L p o = 


1 l 


0 1 


|p o = p o + * 


0 


Because L. generates a line so does T. If a+d = -2 and c / 0 then 
a+i -1 \ / -1 li ja+i -1 1 ^ 

T = 

0; 


a+1 -1^ 

1 1' 1 

M 

/ 

l c °l 

1 ( 0 


1 ( 


The matrix 
If c = 0 then 


/-I 1 

l o -i 


generates a line, or two lines, and so does T. 




T = 


■1 b’ 

0 -1 

which generates one or two lines. 

Case (c) : If | a+d | > 2, consider X, an eigenvalue of T, 


\ a+d , . 1, . ,,2 i9 

X = -y +i \J 1 - (a+d) = e 
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where cos 9 = (a+d) and 

= V*"-I 


sin 0 


(a+d)' 


The sign of sin0 is chosen to be like the sign of b, because of a 
reason that becomes clear later. Note that b / 0 because b = 0 implies 
ad = 1, which implies |a+dj > 2. We will show later the existence of 
a real matrix Q such that 


T = QT C Q = Q 


f cos 0 sin 0 
sin 0 cos 0 


Q 


As generates the circle C, T generates a linear transformation of 
it, which is the ellipse QC. We will construct the matrix Q. As 
and T do not determine Q uniquely, we can expect some freedom in 
the solution for Q. 


QT Q 

c 


-1 


X y \ 

f 

\z w j 

1 l 


cos 0 sin 0 ^ 
-sin 0 cos 0 i 



= T 


Equate components: 

(Ej): a - (xw - yz)cos 0 - (xz + yw)sin0 
(E^): b = (x + y)sin0 
(E^): c = -(z + w)sin0 

(E^): d = (xw - yz)cos 0 + (xz + yw)sin 0 
(E^): 1 = xw - yz . 

2 2 2 2 2 2 

The identity (xw - yz) + (xz + yw) = (x + y )(z + w ) shows that 

one of the equations E-^ through E^ can be discarded, and we can im¬ 
pose another external condition on the system. We choose y = 0, then 


I b 

V sin 0 


w = 


sin 0 Vb sin 0 

L - sin 0 
X Vb sin 0 
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z = 


d-a 


d-a 


2x sin 0 


2V 


b sin 8 


Note that b sin 8 > 0. We get: 

/ 2b 


Q = 


1 


V 


b sin 0 


° 

d-a 2 sin 0 


This is the real matrix Q such that T = QT Q , hence 


2b 0 

; d-a 2 sin 0 


c 
2b 


cos 0 sin 0 
^-sin0 cos 0/ \d-a 


0 

2 sin 0 


-1 


This completes the proof of the theorem. 


(III. 6) On the Spacing Between Successive Points 

2 2 

Theorem : Let us consider the circle x + y = 1, and the hyperbola 

2 2 / 1 \ 

x - y = 1, both passing through the point = [^1 . The 

circle is generated by the matrix T^,, and the hyperbola by 

T„ where 
ri 


T c = 


COS 0 
sin 0 



, and Tj_j - 




The point P on the circle is 
n 


P n = T C nP 0 = 


f cos(n0) 

^ sin(n0) 

and the point P^ on the hyperbola is 
f ch(n/) ^ 

^shfnjzO i 


P = T n P = 
n H 0 


Let d be the distance between P and P , , . 
n n n+1 


We want to 

show that d^ is constant for the circle, but is an increasing 
function of |n| for the hyperbola. 
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Proof: For the circle: 



d = (x , 
n n+1 


- x / + <V„ + 1 - »/ 


Z Z 

= [cos(n+l)0 - cosn6] + [sin(n+l)0 - sin n0] 

r 0 . 2n+l „ . 0-, , r-i 2n+l „ . 0 n Z 

= [ - Z sin— 2 — 0 sin-^ J r [2 cos — ^“ sin J 

, . Z 0 r . Z Znhl „ , Z Znhl ,.20 

= 4 sin y [sin — x — 0 r cos — x — 0J = 4 sm y 


which is independent of n. 


QED 


For the hyperbola P = 


n 

2 , v 2 


;h(n^)^ 

sh(n^) i 


d „' = < X n+l 


x -> + ( y n +i ■ y-> 


n 


n 


= [ch(n+l)^ - chn/]^ + [sh(n+l)^ - sh 


= [2 sh 
= 4 sh 


2r ^ T-- sh ^] 2 + [2 ch 2r ^ + ^ 4 sh 


^]2 
2 


2 ^[sh 2 ^ + ch 2 


2 1 

= 4 sh 2 ^ ch(2n+l)^ 


which is an increasing function of jn| . 


QED 


This theorem is illustrated in Figures III. 10. 1 through III. 10.4. 


Theorem: The Areas Law: Let P = T P^ and let T = det(T). Let S. 
- n 0 l 

be the triangle whose vertices are the origin and the points 


P and P .,, 
n n+1 


Let A be the area of S . Then, A = T n A n . 


n 


O’ 


Proof: 


A 0 ~ 2^ x 0 y l " x l y 0 } 2^ x 0 y 0^ 



=2 P 0 GP 1 
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1 * 

Substitute = TP^ and get = -^PqGTPq. Similarly 

1 ❖ n * n 

A n = -| P o (T ) (GT)T n P Q . 

Consider 


t'gt 


a 

<0 

| 

\ b 

d i 

1 





-ac + ac 

-be + ad \ j 

f o 

"I 

\-ad + be 

-bd + bd j 1 


°/ 


tG 


and (T ) n GT n = r n G. Substitute above and get 

A n ■ lP^(T n fGT"]TP 0 = i T " P ; GTP 0 = t"A q 


QED 


The reader should notice that although there is a similarity this 
is not the Kepler area law. 

Corollary : On a conic all the areas {A^} are equal. 

Corollary : The distance between successive points on ellipses gets 
smaller when the distance of the points from the origin gets longer. 

It is clear that by stretching a uniformly spaced circle into an 
ellipse the uniform spacing is changing into a good spacing. 

The first theorem shows that it is possible to generate conics in 
good spacing, but because of the second theorem any matrix which 
generates a conic, must generate it in good spacing. 

(III. 7) Straight Lines as Degenerate Conics 

There are two kinds of straight lines which are degenerate conics. 
The ones which pass through the origin, and the ones which do not. 

The lines which pass through the origin are asymptotes to hyper¬ 
bolas, and are generated by the same matrices which generate the 
hyperbolas when the initial points are eigenvectors of the matrices. 
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The matrices T, which generate hyperbolas satisfy det(T) = 1 
and trace(T) > 2. These conditions guarantee the existance of two 
distinct eigenvectors, which introduce the two asymptotes. 


The lines which do not pass through the origin, are arcs of an 

ellipse whose major axis is infinite. For example, the ellipse whose 

axes are of length infinity and of length p, is the two parallel lines 

2 2 

y = i\x. The implicit form of the ellipse is y = p . Using the method 
of (III. 3), and substituting a=b = 0- c = 1 (e = -1) we get: 


T(k) - 


1+k 


h* k 2 

*). 

h 

: ) 

0 

1+k J 




As expected, trace[T(k)] = 2. 


This matrix, with the initial point = (xq, generates the 

line y = y^, in the positive direction along the ellipse, which is to 
the right for y n > 0 and to the left for y < 0. All the points of the 
form Pq = (Xq, 0) are eigenvectors of T, corresponding to the eigen¬ 
value 1, and therefore are not changed by the iteration. It is easy to 
show that T does not have any other eigenvectors. 


The generating matrices of proper ellipses [det(T) = 1, |trace(T)| < 2] 
do not have any real eigenvectors, and cannot generate any lines. 


(III. 8) On Curves Created by Matrices with a Non-unit Determinant 

This section suggests, implicitly, building special hardware for 
the iteration P(n+1) = T P(n) or AP(nll) = T AP(n). This hardware can, 
as shown before, generate conics. However it is interesting to know 
what happens if one iterates with a matrix with non-unit determinant 
(which is a necessary condition for conics). In this subsection we 
answer this question. 


It should be noted that if T has negative eigenvalues then the 

g 

function T is not well defined and required an extra care for getting 
continuity. However by observing only integer powers of T, most of 
this danger is bypassed. 
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(III. 8. 1) Consider the cases r = det(T) > 0. The matrix S - t T 
satisfies det(S) = 1, and S generates some conic. If S generates an 
ellipse then T generates an elliptic spiral which winds outward to 
infinity if r > 1, or winds inward to zero if r < 1. See Figures III. 10. 5 
through III. 10. 8. 


If S generates an hyperbola, then T has two distinct eigenvalues, 
X and p. If both are positive then we can define a by \ a = p. Then 


T = H 


H 


-1 


which shows that T generates a linear transformation of y = kx a . If 
both are negative then -T has two positive eigenvalues. T generates 
points which alternate between the curve which is generated by -T 
from Pq, and the curve which is generated by -T from -P^. Note that 
a = -1 implies an hyperbola, and a = 0 (i. e. , 1 is an eigenvalue of T) 
implies a straight line. See Figures III. 10. 9 through III. 10. 11. 


If S generates a straight line then if t / 1, T generates lines 
which belong to a family of curves, all of which pass through the ori¬ 
gin, and go to infinity. The shape of these curves is illustrated in 
Figure III. 10. 12. T generates these curves (or a linear transformation 
of them). Note that one straight line belongs to this family. 


(III. 8. 2) Consider r - det(T) = 0. If T is a non-zero matrix 
then dim[Range(T)] = 1, which means that all {P(s)} belongs to a 
1-dimensional space. Hence T generates a straight line through the 
origin. 

(III. 8.3) Consider t = det(T) < 0. If det(T) < 0 then det(T 2 ) > 0. 

2 

Hence T generates one of the curves discussed before. Every second 

point {P(2n)} which is generated by T is on the curve which is generated 

2 

by T . The other points {P(2n+1)} are on the curve which is generated 

2 

by T through P(l). Hence T generates a sequence of points which 
oscillate between two curves of the same type. 
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(III. 9) Programming Aspects of Coding the Linear Difference Scheme 


This section suggests an incremental method for generating curves. 
We pointed out in (III. 2) that the scheme which is used for generating 
the {P(n)} can be used for generating the (AP(n)}. 

This iteration can be implemented either by conventional pro¬ 
gramming or by hardware. We give in this subsection some "coding- 
tips 11 for programming the linear differences scheme. 

(III. 9. 1) When P(n+1) = T P(n) is coded there is no need to store 
the array of {P(n)} 7 as one can do only with the current value of P ? 
i. e. , the values of x and y. The straightforward coding of the iteration 
is: 


ax + by -> temp 

cx + dy -> y 

temp -> x 


The need for the temporary storage "temp" rises because x cannot be 
changed before it is used for the y calculation. However, the iteration 
can be defined such that x(n+l) is expressed by means of x(n) and y(n), 
but y(nf-l) is expressed by x(n+l) and y(n). 


Consider the following identity: 


fa b ^ 

I 

l c d j 

\ 


Multiply by P{n): 


1 

c/a 

a b\ 
c d j 

1 

^c/a 

This may be coded as: 


0 


(ad-bc)/a 


a b 

0 1 



n 


V y n 


(ad-bc)/c 



ax + by 
ax + Py 


• x 

y 
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where a = — and 0 = ac ^ . If a = 0 but d / 0, a similar scheme works, 
a a ’ 

This eliminates the need for the temporary storing. 

(III. 9. 2) The well-known iteration: 


y 


- 6. 

+ 6 - 


y 


X 

y 


generates an ellipse, because y uses the "new" x as set by the first 
statement. This iteration can be formulated as: 

x , i = x - 6. y 
n+1 n 7 n 


y n+l =y n +e ' x n+l = (l ' 6 > y n + 6x n 


or 


F(n-t-l) = 


( y /n+l 


-6 


1 - 6 * 


= E P(n) 


det(E) = 1, trace(E) = 2-6 <2, hence E generates an ellipse. 

(III. 9. 3) The iteration: 

x -h 6- y -> x 

y + 6* x -> y 

generates a hyperbola, because it is equivalent to: 

1 6 


P(n+1) = 


l+6‘ 


x 


= HP(n) 


\ 


In 


and det(H) = 1, trace(H) = 2 + 6 > 2. 

(III. 10) Examples 

All the pictures (except III. 10. 2) which are appended to this sub¬ 
section were taken from a PDP-1 program which generates curves 
using the method which is discussed in this section. 

Figure III. 10. 1: A circle generated by 16 segments. Note the uniform 
spacing. 
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Figure III. 10. 

Figure III. 10. 

Figure III. 10. 

Figure III. 10. 

Figure III. 1 0. 

Figure III. 10. 
Figure III. 1 0. 

Figure III. 1 0. 

Figure III. 10. 

Figure III. 1 0. 

Figure III. 10. 


2: A circle which was generated according to the method 
which is discussed in Section II. Note the non-uniform 
spacing. 

3: 2 ellipses. The generated points are marked by 

asterisks. Note the n good M spacing, regardless of 
the starting point. 

4: A family of hyperbolas, all of which were generated 
by the same matrix. Note the good spacing. 

5: A circular spiral which is generated by a circle¬ 
generating matrix, multiplied by a scalar, which 
has a magnitude less than 1. 

6: An elliptic spiral which is generated by an ellipse- 
generating matrix, multiplied by a scalar whose 
magnitude is less than 1. 

7. A star is generated by a rotation matrix, with 0 = 277/5. 

8: A star spiral is generated by the matrix of picture 7, 
multiplied by a scalar whose magnitude is less than l. 

2 

9: The parabola y = x is generated by the matrix 

2 

T ~ diag{a,a }. Each part of the parabola has to 
be generated separately. 

3 

10: The cubic y = x is generated by the matrix 

T = diag{a,a } . Each part of the cubic has to be 
generated separately. 

11: The sequence of points 1-2-3-4-5-6-7-8-9-10 was 

generated by the matrix -H, where H is the generating 
matrix of the hyperbolas 1-3-5-7-9 and 2-4-6-8-10. 

12: A family of curves which is generated by a matrix 
which has two equal eigenvalues, but only one inde¬ 
pendent eigenvector. As discussed before, this 
family contains one straight line, which is in this 
example the X-axis. 


43 
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Figure III.10.10 









Figure III.10.11 
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SECTION IV 


INCREMENTAL METHODS FOR HIDDEN-LINE ELIMINATION 

(IV. 1) Introduction 

Producing pictures with hidden lines elimination (HLE) in real 
time is one of the biggest challenges in computer graphics. For some 
years there have existed some programs for generating pictures with 
HLE, like [7], [8], [9], [10], [11], [12], [13], and [14]. All of them 
are many orders of magnitude away from real time computation. 

There exists only one system which produces images, with HLE in 
real time. This is a special-purpose system which was designed and 
built by GE for NASA, [15], at a cost of about $3, 000, 000. 

Recently, John Warnock of the University of Utah devised an 
algorithm [16], which for the first time brings some hope for econo¬ 
mical real time HLE. The programs mentioned before use a 
straightforward brute force algorithm which checks each possibly-seen 
entity against each possibly-hiding entity. This checking of "all 
against all 11 makes the required amount of computation to be propor¬ 
tional to the square of the number of defined objects. The Warnock 
algorithm (WA) deals with the objects according to the order in which 
they are located in the picture, not according to their arbitrary order 
in the data structure. The amount of computation required by the WA 
grows at a rate less then the square of the complexity. In (IV. 2) we 
describe briefly the WA. In (IV. 3) we show an incremental approach 
to the WA which eliminates redundant computation by organizing the 
computation in such a way which saves at any step the information 
which can be used in later steps. It is estimated^ that this approach 
can cut the required amount of computation for a typical simple figure 
(like the picture of a house), by an order of magnitude. A most time- 
consuming problem, which is at the core of most HLE programs, is 
finding whether a given point is inside of a given polygon. In (IV. 4) 
we show an incremental method for solving this problem. 


1 


By Mr. 


Warnock and others. 
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(IV. 2) A Brief Discussion of the Warnock Algorithm 


The WA has two basic logical units, a "control-unit" and a 
"looking-unit n . The control unit chooses a portion of the picture, 
which I call a window, and tells the looking-unit to work on it. 

The looking unit considers this specified subpicture (the "window") 
and finds what is seen in this window, or announces a failure to find 
it, due to a too complex situation. In case of success, the control 
unit outputs the results to some display system. In case of failure, 
the window is put on a list of "unsolved-windows". Later each un¬ 
solved window is subdivided into some smaller windows, each of 
which is given to the looking-unit for consideration. The looking-unit 
never announces failure when the size of the window has been reduced 
to a single resolution unit of the display. This guarantees that the 
algorithm is always completed in a finite number of steps. The pro¬ 
cess continues until all the portions of the picture are solved. 

There are many possible versions of the algorithm, depending 
on the complexity of the looking-unit and the control-unit. The simplest 
looking-unit is one which announces failure if any vertex (point) or 
edge (line) is seen in the window, i. e. , if its projection lies inside 
the window, and no opaque surface is between it and the observation 
point. If there is some surface which hides everything else in the 
given window, then the looking-unit tells the display system that the 
"color" of the window is the "color" of this polygon. If nothing is seen 
in the window, then the looking-unit announces the window to have the 
"color" of the background (mostly black). The simplest control-unit 
is the one which always divides windows into four identical quarters. 

A more sophisticated looking-unit can handle more complex 
situations, like a window with a single line in it. Handling these win¬ 
dows, without announcing failure saves a considerable amount of com¬ 
putation; the more sophisticated the looking-unit is, i.e., the more 
complex situations can be handled without announcing failure, the less 
subdivision is required. 

The simplest control-unit always subdivides windows in the middle 
however, it is advantageous to subdivide windows, at some complex 
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point, like a vertex. By subdivision at a complex point, the com¬ 
plexity is most often reduced to a simpler case. 

There is a wide range of looking-units and control-units which 
can be used with the algorithm. A beautiful property of the algorithm 
is the independence of the structure of the control-unit and the looking- 
unit . 

(IV. 3) The Incremental Approach to the Warnock Algorithm 

The basic idea of the incremental approach is reordering the 
computational procedure in such a way that essential information can 
be saved, and used in later steps. Storing this information eliminates 
the need to ask the most time-consuming questions more than once. 

Let W = (Wp W^, W^, W^} be a window which is divided into 
the 4 windows: W^, W^, W^, as shown below: 


W 1 

w z 

W 3 

w 4 


Figure IV. 3.1 : A window subdivided into 4 subwindows 

Let the index 1 always denote the upper left corner, 2 for the upper 

right corner, 3 for the lower left and 4 for the lower right one. If 

W. is subdivided we have W. ~ {W., W. 9 . W..}. If W. . is 

i i il’ i^ ? i3’ i4 ij 

subdivided we have W . = {W. .. I k = I, 2 3,4} and so on. 

ij ijk' ’ ’ ’ 
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Let us consider the following window, W: 


W 1 

W 21 





W 23 

W 

24 

w 

31 

W 32 

W 4 

33 

W 34 


Figure IV. 3. 2: Subwindows divided into subwindows 


The window W, with its subdivisions, as shown in the above figure, 
can be represented by: 

W = {W lf W 2 , W 3 , W 4 } 

W 2 = ^ W 21 ’ W 22> W 23’ W 24^ 

W 3 = {W 31’ W 32> W 33’ W 34 } 

W 22 = ^ W 221 ’ W 222’ W 223’ W 224^ 

By substitution we get the following representation: 

w = {w p {w 2 l ,{w 221 ,w 222 , w 223> w 224 },w 23 , w 24 }, 

{ W 3P W 32 , W 33 , W 34 ) , W 4 } 

The graphical representation of the above is as follows: 


59 














Figure IV. 3. 3 : A tree representation of Figure IV. 3. 2 


Each node {W.. } corresponds to a window. A window which is 

ij. . . 

declared to be too complicated is represented by a node with 4 sons, 
a solved window is represented by a terminal node. The task of 
the control-unit is to decide which node should be considered next. 
Tree scanning can be done in many ways. One of these is the 
"constant-depth walk M *. In this example, a constant-depth-walk 
visits the nodes in the following order: 


{w, w,, W- 


w w 

3’ 4’ 


W 


W 


W 


W 


W. 


W, 


w. 


21’ 22’ 23’ 24’ 31’ 32’ 33 


W 


W 


W 


W- 


34’ 221’ 222’ 223’ 


W 224> 


This order is achieved by starting with W, and keeping a FIFO (First 
In, First Out) stack for the unsolved windows. This ordering implies 
that windows are processed in the order of their sizes. 

The incremental approach is to scan the problems-tree by a 
2 

prefix-walk , which is achieved by keeping the unsolved problems in 


1 See [17]. 
^See [17]. 
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a LIFO (Last In, First Out) stack, instead of a FIFO stack. In our 
example, the prefix walk is: 


{w wwwww w w w ww 

\w, W z , w 21 , w 22 , vv 221 , vv 222 , vv 223 , vv 224 , w 23 , w 24 , 


W 3> W 31> W 32> W 33’ W 34> W 4 } 


The length of this LIFO stack cannot exceed N, the maximal number 
of subdivisions which is the base-2 logarithm of the scope size (e. g. , 
N = 10 for a 1024x1024 scope). There are never more than N open, 
unsolved problems, which reduces the amount of storage required. 

The most time-consuming problem is finding the relation 
between a given polygon and a window. This relation is one of the 
following: 

(a) the polygon is outside the window, 

(b) the polygon surrounds the window, 

(c) some edges of the polygon intersect the window, but no 
vertex is inside the window, 

(d) some vertices of the polygon are inside the window. 

If the polygon P is outside the window W (property (a)) then P is out¬ 
side any subwindow of W. 

If P surrounds W, it also surrounds any subwindow of W. 

If P has no vertex inside W, it has no vertex in any W., but P may 

be outside some W., or surrounds some W.. 

i* J 

If P has any vertex inside W, then all the above relations can hold 

between P and W.. 

i 

This can be represented in the following diagram: 
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Figure IV. 3. 4 : The relations between the polygon/window relations 


Outsideness and surroundedness are terminal because they do not 
change for descendants. When any terminal-relation exists between 
a window W and a polygon P, this relation is known to hold between 
P and any subwindow of W. The classification of the relations into 
terminal and non-terminal is the key to the incremental approach, 
and to the prefix-walk. When a relation is found between a window 
W and a polygon P, the relation is stored. If this is a terminal rela¬ 
tion, then it is known to hold between P and any descendant of W (i. e. , 
subwindow). This saves the repeated checking for the relation be¬ 
tween polygons and windows. If, when a polygon is considered, there 
is no terminal relation between it and any ancestor of the window, 
then the program should look for it. The relation between windows and 
polygons are stored in a stack which is associated with the polygons. 
The length of the stack is N (as defined before) and each entry can be 
expressed by 2 bits. Hence, by associating a stack of size 2N bits 
with each polygon we can convey its relation with any window, to all 
of its descendants. This works only for prefix walk, and not for 
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constant-depth walk, because the number of open-problems on constant- 
depth walk can be as high as 4 , and associating 2- 4^ bits to each 

polygon is very likely to be infeasible. The prefix walk contains never 
more than N open-windows. A flow chart of the algorithm, using the 
simplest looking-unit and the simplest control-unit is shown in 
Figure IV. 3. 5. 

(IV. 4) The Incremental Solution for the In/Out Problem for Polygons 

In the core of all HLE programs there is a solution for the in/out 
problem. The computation required for solving this problem is a great 
portion of the total computation required for HLE. It is very important, 
therefore, to improve the technique for solving this problem. We will 
show an incremental solution for this problem, and outline the hard¬ 
ware implementation of this technique. 

The In/Out problem for polygons can be stated by: 

(Ql): "Is the point P^, inside the polygon P?"^ 

Another form of this problem is: 

(Q2): "Does the polygon P surround the point Pq?" 

A generalization of (Q2) is the following: 

(Q3): "How many times does the polygon P surround the point P^? n 

If the edges of the polygon P do not intersect themselves, then 
(Q2) and (Q3) are equivalent. But this is not necessarily the case. 

There are two techniques for solving this problem. One is based 
on angles summation and the other on intersections of the polygon and 
a line from the point P^ to infinity. The angle summation technique 
works as follows. Define: 0. = £ P.P P . where -7T < 0. ^ 7T . Compute 


An N-vertices polygon is given by the ordered set of its vertices 
{P^ | i = 1,2... N} . We define P^_^ = P^, suc ^ that the edges are 

E. = PTP. , , for i = 1, 2. . .N. 
i i l+l ’ 
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l‘hc simplest Warnock algorithm for line drawing. 


Figure IV.i. j : 






















N 

CT = 2 0. . 

i=l 1 


Here or is the total angular change of the polygon P around the point P^. 
If P surrounds the point n times then cr = 27Tn. Note that n and cr may 
be positive or negative, depending on the sense of the polygon. Note 
that according to Euler's theorem 


_ _cr__ 

~ 27T “ 27Ti 



Pq is outside P if n = 0, and inside otherwise. Hence or = 0 if Pq is 
inside, and [ cr j > Ztt if Pq is outside. 

In implementing this technique one computes 
N 

or = 2 e. , 

i=l 1 


and announces insideness if ]cr| ^ 7T* and outsideness otherwise. The 

computation of the angles {0.}, is very expensive timewise although 

1 7T 

each 0^ has to be computed only with accuracy of c = ^ • 

The other technique, the ray-method works as follows: we define 
a ray from Pq to infinity. For simplicity in this discussion we assume 
that Pq = (0, 0). Let the ray be R = {0 < x < oo, y = 0} . We check all 
the edges of P to find if they intersect R. Let n be the count of inter¬ 
sections. P^ is announced to be outside, if n is even, and inside 
otherwise. 

The angles - summation method can answer the problem (Q3). The 
ray method can answer only (Q2) for certain cases, but is much faster 
since intersecting the {E^} with R is simpler than finding the { 0^} . 

Consider, for example, a polygon which describes the digit "4", 
as illustrated in Figure IV. 4. 1. Consider the points A, B, C and D. 
Applying both methods to these points shows that A is inside but B and 
D are outside. However the two methods do not agree about the point 
C. According to the angle-summation method C is inside (twice) but it 
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Figure IV. 4. 1 : A polygon describing the digit ,T 4 n 


is outside according to the ray method. The technique described below 
is as powerful as the angle-summation method, and is at least as fast 
as the ray method. 

The order of the vertices in the definition of P implies a direc¬ 
tion (clockwise or counterclockwise) for the polygon. This direction 
is used to associate a value (plus or minus one) to each intersection 
of edges and the ray R. The algebraic sum of the values of the inter¬ 
sections is the number of times in which P surrounds the point Pq. 

For simplicity again let P^ = (0, 0). We choose the ray to be 
R = {0 < x < oo, y = 0} . The edge intersects R only if y^ and 
have different signs . The value of the intersection is defined to be 
+ 1 if y. ^ 0, and -1 if y. < 0. If E. does not intersect R the value of 
their intersection is defined to be zero. 

By considering only the signs of (x^, y^) and ( x j|pY-||) we can 
obtain the value of the intersections for the following cases: 


Note that zero is always considered as a positive number. 
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Case 

(a) 

(b) 

(c) 

(d) 

(e) 

(f) 

(g) 

(h) 

(i) 


x. 


$ 

$ 

+ 

+ 

+ 


4 > 

+ 

+ 

+ 

+ 

+ 


i+1 


3> 

$ 

+ 

+ 

+ 

+ 


y i+l 


4> 

+ 


+ 

+ 


V(V = value) 

0 

0 

0 

+ 1 
-1 

+ 1 or 0 
+ 1 or 0 
-1 or 0 
-1 or 0 


Cases (a)-(e) cover 12 out of the 16 possible sign combinations. The 
other 4 combinations are (f)-(i). The intersection values in these cases 
have to be found by using arithmetic operation. A possible set of rules 
which resolves these cases is the following. 

Cases (f) and (g): if ^^2 < x 2 y l t * len V = +1, otherwise V = 0. 

Cases (h) and (i): if x^y^ > x 2^2 ^ en ^ = "1, otherwise V = 0. 

For the general case, where Pq is not the origin, these cases are 
described by the following Karnaugh map: 

y B = 1 


= 1 


0 

0 

? 

+1 

\ 

0 

0 

0 

? 

\ 

? 

0 

0 

1 

0 

-1 

? 

0 

/ 

0 


X A = 


= 1 


X A = 1 

if 

X. < 

1 

x o 

II 

< 

if 

y i < 

y o 

X B = 1 

if 

x i+l 

O 

X 

V 

y B = 1 

if 

y i+l 

A 

o 


X B = 1 


Figure IV. 4. 2 : A Karnaugh map for the signs conditions 


1 


for don't care. 
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The (?) indicate the cases in which more computation is required for 
finding V. 

The In/Out problem as described before deals with the relations 
between a polygon and a point. The WA deals with relations between 
polygons and windows. Since a point may be considered as a degenerate 
window, solving the relations between windows and polygons is a more 
general problem then solving the relations between polygons and points, 
as described before. The technique in this case is to check all the 
polygon vertices {P^} for intersections of the edges {E_.} and the window 
sides, and for being inside the window, W. If some vertices are inside 
W then this is a case of complexity (d) as defined in II. 3. If no ver¬ 
tices are inside W, but some edges intersect W then this is a case of 
complexity (c). If no vertex is inside W, and no edge intersects W then 
this is either case (a) of the window being outside the polygon, or case 
(b) of the polygon surrounding the window W. It is still to be resolved 
whether this is case (a) or (b). 

The same relation (outsideness or surroundedness) which holds 
between the polygon and the window W holds also between the polygon 
and each of the window corners. Hence, it can be found for W by 
finding it for any one corner using the method described earlier. How¬ 
ever, this computation which scans all the {E^} can be combined with 
checking the {E^} for intersecting the window sides, as both computations 
require common calculations and comparisons. 

It is feasible to build a special-purpose hardware, using this 
technique, which scans once all the {p^} and announces which portions 
of the {E.} are inside the window W, if any. If none, it announces an 
outsideness or surroundedness relation. This n box n can get the re¬ 
quired information {P^} by using a direct channel to some memory 
and work in parallel and independent of the host computer. 

This hardware is primarily a clipping-divider [19] which clears 
some n-register when it starts working on a new polygon, and updates 
it when working on the polygon edges. After the last edge has been pro¬ 
cessed, the n-register contains the number of times in which the polygon 
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surrounds some window corner. If there are no intersections then 
the n-register shows how many times P surrounds W. For each edge, 
E^, the clipping divider compares the coordinates of and p m with 
the window coordinates. These comparisons are the inputs needed 
for updating the n-register according to the table of Figure IV. 4. 2. 

In the 4 complicated cases (which are indicated by - ?- in the table) 
the clipping-divider cannot make a trivial decision either and finds a 
point P^, on the edge E^ ? which together with each vertex is used for 
two simple decisions. This point P^, enables the In/Out logic to make 
two decisions according to the table. This may be illustrated by the 
example in Figure IV. 4. 3. 



Figure IV. 4. 3 : Examples for window/line relations 
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SECTION V 


PRODUCTION OF HALF-TONE IMAGES IN REAL TIME 

(V. 1) Introduction 

Most computer display systems use calligraphic (line drawing) 
scopes. Real world pictures (like the ones produced by cameras and 
TV) are half-tone images. Production of real time half-tone images, 
at a feasible cost'*' has been for a long time, a challenge for computer 
graphics. 

The reason why this challenge was not met for a long time, is 

2 

the huge number of points which have to be intensified for displaying 
a half-tone image. The bigger this number is, the less time is available 
per point. The short time which is available per point, might be enough 
for intensity calculation, but not for random beam positioning. For this 
reason, any half-tone image production must be based on some standard 
scanning pattern with the beam intensity-modulated according to the 
image to be displayed. 

Using a scan technique for displaying introduces a sorting problem. 
The scanning requires a certain order in which the data must be arranged, 
while the data is generated by the computer according to some order 
which does not necessarily agree with the scan order. This need for 
sorting creates what is known as the "scan-conversion 11 problem. In 
this section we describe two methods for producing a half-tone image. 

One, which is based on a TV scan by using Cartesian coordinates, is 
described in (V. 3). The other method, which is based on a radar scan 
using polar coordinates, is described in (V. 4). 


'''There exists a system in the MSC (NASA) in Houston, Texas, which 
was built by GE for about $3, 000, 000. [15]. 

2 

For a 512x512 raster scope, with refreshing rate of 30 times per second, 
about 7. 5* 10 6 points have to be intensified per second, which is about 
133 vsec per point. 
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Both techniques use the same two functional units, a line-controller 
(LC) and a line-controller-interface (LCI). The LC generates the 
scanning and modulates the beam intensity according to its position. 

The LCI prepares information for the LC. There is a trade-off between 
the complexity of the LCI and program which feeds it. The simpler the 
LCI is, the more has to be computed by the program, and vice versa. 

The main feature which makes these techniques feasible is supplying 
information only about the points in which the intensity rule changes 
along scan lines, instead of supplying information about each point. 

In (V. 2) below we describe the general LC, which might work for 
any 2D-coordinates system. 

(V. 2) Line Controller for Generalized Coordinates 

As the scope face is a 2D surface, there are many ways to map 
it by using some 2D-coordinate system, which is not necessarily car¬ 
tesian. 


Let (u, v) be the coordinate system which is used to cover the 
scope face. Without loss of generality we can assume that the displayed 
area is: 


{(u, v), 0 < u < 1 3 


0 « v < 1} 


This coordinate system is useful only if there exists a fast technique 

for generating the scan lines: {(u, v.), 0 ^ u ^ 1} for all values of v.. 

J J 

Let us assume that this is the case. 


A (u, v)-scan image is the following sequence of intensities: 
{{i(u, v) | 0 ^ u < 1} , 0 < v < l} ^ 

The image is generated by generating the sequence 
{v 0 = 0, v p v 2 v N = 1} 

synchronized with displaying a frame, and supplying the scan lines: 


If the scanning is generated in an analog technique, the LC is syn¬ 
chronized with the scanning. 

2 

i hereafter is for intensity 1 . 
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{i(u, v n ), 0 < u < 1} 

for each v . 

n 

t 

The LCI prepares a list of positions at which the intensity rules 
change along the scan line v = v . For simplicity we assyme that all 
intensity rules are just constant intensity in some area. We will des¬ 
cribe other intensity rules below. An intensity profile like the following: 


I 



Figure V. 2. 1 : Intensity profile 


is described by the following list: 

(0, i B )(u l5 i c )(u 2 , i A )(u 3 , i D )(u 4 , i A ) . 

* 

The intensity which is associated with each m is displayed between 

u = u. and u = u. - (if there is no u. , , , then u. , , = 1). 

J J+l J+ 1 J+ 1 

When the LC starts displaying a frame it clears its v-register (VR) 
its u-register (UR) and loads its intensity register (IR) from the first 
piece of data. The UR and the VR are used for controlling the beam 
position and generating the scan, the IR modulates the displayed intensity. 
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The next data (u^,T) is loaded into the next-u-register (NUR) and into 
the next-intensity“register (NIR). For displaying each point, the UR 
is incremented by some 6U, and the IR sent to the scope. If the in- 
cremented UR is equal to the content of the NUR, then the NIR is 
jammed into the IR for the next point, and new data is entered into 
the NUR, and the NIR. When the UR reaches the value u = 1, the 
VR is incremented by 6V, and the UR and the IR are loaded from the 
NUR (which contains u = 0) and the NIR, while new data is placed in 
them. When the VR reaches V = 1 it is the end of the frame. Note 
that the magnitude of Su determines the resolution along scan lines, 
as well as the scan velocity, and the magnitude of Sv determines the 
number of lines per frame. 


The organization of this LC is illustrated in Figure V. 2. 2. If 

the shade between successive data points {u.} is not constant, but 

0i ^ 

changes linearly with u, (i. e. , -r— is constant) then 2 more registers 
0i cfu 0^ 

are needed, the - register (JUR) and the next - - register (NIUR). 

The organization of this L*C is illustrated in Figure V. 2. 3. 


The generalization for higher degrees of shading rules is obvious, 
and would follow the polynomial generation scheme described in (II. 4). 

Note that if 6U is of the form 2 then adding 6U is only a 
counting operation, and the U = 1 condition is recognized by a carry 
from the most significant bit, or by the change of the UR to zero. 

When this condition is detected, there is no need to clear the UR by 
loading the NUR to it, as the binary representation of u = 1 can be 
u = 0. The above remarks hold also for 6V. 


Note also that the process of scanning which involves repetitive 
addition of the SU to the UR, and of the SV to the VR is a simple inte¬ 
gration process which can be carried out by analog means. 

(V. 3) Production of TV-scan Images 

The technique described below was first suggested by Professor 
D. Evans of the University of Utah. The version which is described 
below, and which was simulated on a PDP-1 computer is a result of 
many discussions with Professor I. Sutherland. 
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I'lyuro V.2.2 : A shader for constants intensities 
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The TV-scan is based on cartesian coordinates where the (u, v) 
of (IV. 2) are replaced by (x, y). The LC is exactly as described 
there, however, the LCI can take advantage of the properties of 
cartesian coordinates, and perform some operations which otherwise 
have to be done by the program. 

The LCI is able to accept a description of a linear shading rule 
line, and produce by itself the data about all the intersections of this 
line with scan lines. This is done by the following scheme; consider 
the line which initiates^" the following linear shading rule: 

i(x,y) = i(xj,yj) + (x - Xj)|^ + (y - • 

For the scan line y - y^ where y^ ^ y^ ^ the L»C needs the following 
information: 


( 3i . 

x n’ l n’ 


where x = x, + (y - y, ) 


x 2 - x } 


n '1 ‘ M y 2 - yi 


which is the x coordinate of the inter¬ 


section of P,P n with the scan line y = y , and i is: 

1 2 3 3 n y n 

V = i<X n’ y n' = i<x l- y l> + (x „ ' x l’to + <y n ' y l 1 '" 


3y 


which is the intensity at (x^, y^), the intersection point. 

For the next scan line y - y = Y n + ^Y the LC needs the following 
information: 

( 9ix 

' X n+1’ n+1’ 3x 

where 


n+1 


= x i + ( vn ■ 


x 2 " X 1 


1 y? "Yi 


■ * dx 

- X + Oy - — 

n 7 dy 5 


and 


i. e. , this shading rule applies on the right of this line, while the 

scanning goes left to right. 
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n+1 




+ (x. 


n+1 


x 9i 

- x, ) —— 
i' ax 


+ (y 


n+1 


v 8i 

• y i’i7 


i + 
n 


i £1 dx , j)i_ 
9x dy 9y 


)6y = 


i + 4“ 

n dy ; 


Hence, the LCI can generate the {x } and the {i } for successive scan 

dx A' n 

lines if it is given 6y * and 6y • ^ , in addition to the data about P . 

the first intersection with the scan lines. It needs also some infor¬ 
mation which indicates where P ^ is, to terminate this shading rule. 


To summarize: The LCI can feed the LC with all the needed 
data for a shading rule boundary if it has the following information for 
each shading rule boundary: 


(x 


1>?V 


dx 

dy 


cU 

dy 


6 Y,i(x 1 ,y 1 ),||,n) 


where ( x pyj) is the first intersection of the shading boundary with the 
scan lines. 


dx 

dy 


indicates the slope of the boundary 


6y 


dl_ 

dy 


the intensity change along the boundary , as a result 
of changing y 


the intensity at ( x pYj) 

the intensity change along the scan line , as a result 
of changing x 

the number of scan lines which intersect this boundary. 

From this information it can generate for each scanning line the 

(x , i , —) which are needed by the LC. 

n’ n’ 3x 

The usefulness of this LCI depends on the nature of the picture. 

If it has many long lines which are intensity rule boundaries it pays 
off to build such an LCI, otherwise a simple output channel from the 
computer storage can serve as a LCI. 


KXi.y.) 


M. 

9x 


N = 


^2 ' Vl 
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(V. 4) Production of Radar Scan Images 


The radar scan is generated by using the polar coordinates (r, </>) 
over the domain 0 < r ^ 1 , 0 ^ ^ < Ztt . The scan lines are the radii 
{(r,/), 0 < r ^ l}, for all 

2 

We will describe below a LC which uses the (r , /) coordinates, 

but generates the (x, y) coordinates for the scan, as required by the 

2 

orthogonal deflection mechanism of the CRT's. It uses r rather than 
r, because mathematically it is easier to find the square of a distance 
between two given points rather than finding the distance itself. 

The operation of this L»C is along the lines described in (V. 2) 

2 

before, where u is replaced by r , and V by $/Zl7. 

The (x, y) coordinates for (r, $) are: 

x = r cos $ 
y = r sin / . 

The cos / and sin^ define the direction of the scan line, and are 
supplied by the LCI for each scan line. The 2 multiplications which 
are needed for the x and y calculation can be replaced by two additions 
because r is incremented linearly. The (x, y) generation can be gen¬ 
erated by a scheme such as the one shown in Figure V. 4. 1. Another 
possible implementation is shown in Figure V. 4. 2. 

Note that the r-register has no function in this implementation. 
Note also that if the r-register has m bits, and the cos / and sinfi 
have n bits then the x and y registers need m + n bits. 

2 

Next we have to generate r for the comparison of the u-register 

2 

with the next-u-register (see V. 2). The {r } can be generated, of 
course, by a multiplication. But this multiplication can be replaced 
by an addition, by using: 

r -1 

r = 2 (2n + 1) . 

n= 0 
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igure V.4.1 : A scheme for generating X and Y. 
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Figure V.4.2 : A multiplication £ roe scheme tor computing X and Y. 
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Note that adding 2n + 1 is implemented by adding n, shifted left one bit, 
and adding 1 to the least-significant bit. 

Incorporating all the above to one system suggests the implemen¬ 
tation shown in Figure V. 4. 3. 


8 1 



I'igure V. A. 3 : An lr', I sluicling system. 
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SECTION VI 


FAST INTERSECTING OF LINES 

(VI. 1) Introduction 

Another problem in computer graphics is the line intersecting 
problem, i.e. , finding if two given line segments intersect, and if 
so, where. 

The mathematics required for solving this problem directly and 
classically is simple, in principle, and requires only some additions, 
multiplications and a division. However for some applications, the 
time required to do these operations precludes solving this problem 
in real time. 

Our goal is to achieve a method for solving this problem very 
quickly which can be implemented by special-purpose hardware. It 
is impossible to separate the method from its hardware implementation. 

In (VI. 2) below we show and describe the operation of the 
"Sutherland interpolator" which is the core of the line intersecting. 

A description of the Sutherland interpolator can be found in [19]. 

In (VI. 3) we show how the Sutherland interpolators are used for 
intersecting line segments. 

In (VI. 4) we describe a device which can get (from some memory, 
via a channel) a string of contiguous line segments, and finds the inter¬ 
sections of these segments with a given segment. 

(VI. 2) The Sutherland Interpolator 

The basic function of the Sutherland Interpolator (SI) is finding 
the intersection of a given segment with one of the coordinates axes. 

Let the segment be P A P B , where P A = (x A , y A^ and P B = ^ x B’ y B^‘ 
Finding the intersection of this segment with the y-axis is finding the 
value y^, such that the point = (0,y^) is on the segment p ^Pg, as 
shown in Figure VI. 2. 1; 
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Figure VI. 2. 1 : A line intersecting the y-axis 


The solution satisfies: 


y B " y 0 


B 


y A ’ y 0 X A 

It is self-evident that is between P A and P_ or P^ e (P A , P_,) if 

0 A B 0 A’ B 

and only if and x^ have different signs, or x^x^ < 0. 
Eliminating y^ from the above condition yields: 

X B y A " x A y B 


y n = 


X B " X A 


Note that < 0 implies that x n - x A / 0. 


B 


The basic operation of the SI is finding P the midpoint, of 
the given segment: 

P = (x m’ y m> Where x m = I (x A + X B> and y m = I * y A + y B> 


If x^ = 0 then y^ = y^ and the problem is solved. Otherwise: if x 


m 


m 


m 


has the same sign as x A then the solution, P^, is between P and 

A 7 0 ? m 

P-d : P^ e (P , P-d ). If x has the same sign as x^ then P n c (P , -P * )« 
B 0 m 9 B' m B 0 m ? A' 
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Note that x A x_. < 0 implies that either x x A < 0 or x x_> < 0, but 
A B r m A m B 9 

not both. 

To summarize: after one iteration of finding midpoint we either 

solve the problem (if x^ = 0) or reduce the problem to a similar one 

with the segment (P , P A ) or (P , P_>) instead of (P A , P_ ). Note that 
6 v m J A 7 v m* B v A ? B 

the width of the new segment is half of the width of the original segment. 

By iterating this operation M times, replacing at each step one 
of the end-points by the midpoint, the width of the segment is reduced 
by the factor 2 For N-bit numbers this process cannot last more 

than N steps until the problem is solved, either by hitting P^ or by 
default when the width of the segment which contains P^ is reduced to 
one unit. 


Each step consumes one add-time . No time is required for the 
division by 2, which is achieved by the way the registers are connected, 
ignoring the least significant bit of the adder, and reproducing its most 
significant bit (the sign bit). 


The total duration of this process for N-bits number is bounded 
by N-add-times (plus a little overhead) which is equivalent to one -divide 
time. 


Before describing the logic and the control of the hardware imple 
mentation of the SI, we make some remarks: 


/ 

if x ^ x 2 > 

Define S{x^, y ^, x^, y^} = < if x^ = x^ = 0: 

, otherwise: 


undefined 

undefined 

x 2 Yl - Xl y 2 



From the definition it follows that the solution of 


is 


yj . ~ y x i 

YZ - Y x 2 


.md y r (y ,, y 2 ) 


y = S{x l3 y l 5 x 2 ,y 2 } 
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It is easy to verify that for any a / 0: 

S{ a x 1 ,y 1 , a x 2 , y 2 ) = S{x 15 y 15 x 2 , y 2 ) , 

As illustrated in Figure VI. 2. 2. 

Y 



Figure VI. 2. 2; X scaling does not change the intersection 


The accuracy of the solution depends on the slope of the segment. The 
flatter the segment is, the less sensitive the solution is. Therefore 
it is advised to normalize the x*s, by scaling them up, before loading 
the X-registers. 

It is also easy to verify that: 

S{ -1, 0, b-1, a} = ^ (for b > 1) 

which suggests a way of using the SI for division. This identity can 
be illustrated by the following figure: 
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Figure VI. 2.3 : Using the Sutherland interpolator as a divider 


If b < 1 but 2*b > 1, then: 

S{-1, 0, 2V1, a} = — = 2~ £ J 

2^b b 

and 

2 £ S{-1,0, 2^-1, a> = 2 £ S{-2' £ , 0,b-2" £ ,a} = g . 

The Operation of the Hardware Implementation 

After loading the Xj,yj, x 2> y 2 registers, and verifying that 
x l x 2 ^ the iterations can start. 

Step (a): Xj + X 2 —> 2 x, and y^ + y^ —>2y. Then continue 
with steps (b), (c), and (d) simultaneously. 

Step (b): If 2 x = 0 then y Q = ^Sy, stop the iterations. 

Step (c): If sign (2 x) = signfxj) then: ^Sx -»Xp |sy- >yj. 

Then go to step (a). 

Step (d): If sign(2 x) = signfx^) then: ^2x-^ ^ Y- *^2' 

The organization of this logic is illustrated in Figure VI. 2. 4. 
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Figure VI .2.4 : A schematic 
<1 raw in}; ol the Sutherland 
Interpolator. 
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Remarks 


(a) The division by 2 is achieved by ignoring the LSB (least 
significant bit) of the sums, and reproducing the sign 
bits. In 2‘s complement arithmetic negative numbers do 
not converge to zero, by repeated arithmetic right shift, 
but to -1. Hence in step (b) it is important to check also 
for Sx = -1. 

(b) If N-bits adders are used for N-bits numbers it is necessary 
to add the logic for detecting overflow. However, as the 
least significant sum bit is ignored, and only the least signi¬ 
ficant carry-bit is used, the last adder stage can be replaced 
by an 11 AND 11 gate, freeing one adder stage for taking care of 
the overflow which can be done by treating the numbers as 
(N+l)-bits numbers, with the sign bit reproduced, which do 
not overflow. 

(c) Checking for signs equality between x^ and can be done 

during the iterations, instead of before them. If sign(x^) = sign(x^) 
then 

[sign(Sx) = sign(x^)]* AND- [sign(£ x) = signfx^)] , 
which is an easily detected condition. 

The Control of the SI 

The SI can be in one of the following states: 

(A) waiting for input in the registers; 

(B) iterating; 

(C) waiting for output from the y-adder, to be recorded. 

State (A) is always followed by state (B). 

State (B) is followed by state (C) only in case of success (intersection) 
or by state (A) in case of failure (no intersection). 

State (C) is always followed by state (A). 

This can be represented by the graph shown in Figure VI. 2. 5. 
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Figure VI. 2. 5 : The states of the Sutherland interpolator 


Let us assign 2-bits codes to the states: 

(0, 0) - State (A) 

(0, 1) - State (B) 

(1, 1) - State (C) . 

The remaining code (1, 0) will be used as a transit state between (C) 
and (A) to insure that (B) does not occur accidentally. The control 
logic can be implemented as shown in Figure VI. 2. 6. 

(VI. 3) Lines Intersector 

The SI as described in (VI. 2) is capable of intersecting segments 
with the Y-axis. In this section we will describe how to intersect any 
two segments, given by their end-points. 

For intersecting P^P^ with a line which is parallel to the Y-axis: 
x = Xq, a simple coordinates translation is required. It is easy to 
verify that if = (x^ ? y^) is the intersection of P^P^ and x = Xq then: 

x 0 = S{x A ■ x 0’ y A’ x B " x 0’ y B } • 

Intersections of segments with the X-axis or with lines parallel 
to it can be found by changing the roles of x and y in the SI. Example: 
the intersection ^^g with y = y^ is Pq = ( x q> Yq) where: 

y 0 = S{y A ' y 0’ x A’ y B " y 0’ y B } • 
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Figure VI.2.6 : The control logic for the Sutherland Interpolator. 
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Next we consider the problem of intersecting P A P. R and P^P^, 


A B 

where neither segment is parallel to the coordinate axes. The 
solution Pq = (x n , y n ) satisfies: 


0 > y 0’ 


x, 


x 


B 


x. 


y A 


and 


B 


X 1 " X 0 


x^ 


X, 


r 0 y 2 


Eliminating x^ and y^: 

(x A y B " x B y A )(x 2 " x l> " (x l y 2 " x 2 y l )(x B " X A> 


X 0 


(y 


B 


y A )(x 2 


' x x ) - (y£ ~ y 1 )(x B - x A ) 


_ ( x A y B ' x B y A )(y 2 ' y l } ' (x l y 2 ' x 2 y l>( y B ' y A> 

Y ° " (y B " y A )(x 2 " x l ) " (y 2 " y l )(x B " X A> 

Needless to say, these equations are not the basis of our technique. 

A possible solution, is reducing the problem into a solved one, 
by rotating the plane until the line P^P^ is parallel to the Y-axis, 
then using the SI to find the intersection in the rotated plane, and 
rotating back to get the true intersection point. 

The disadvantage of this technique is twofold. First it requires 
rotation of the 4 given points plus the solution point. Rotating each 
point requires 4 multiplications. Hence, all together 20 multiplica¬ 
tions are needed. However, all the multiplication required for rotating 
the 4 given points can be done in parallel consuming only one multiply- 
time (but 16 different hardware multipliers). Second, for rotating 
points, we need the sine and the cosine of the angle between the line 
P^P^ anc * coorc ^ nate a *es. 

We shall show how to avoid these difficulties by using only 4 to 
6 multiplications, avoiding the back-rotation and avoiding the need for 
sines and cosines. We will do it by the following scheme: 
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Consider the example in Figure VI. 3. 1: 



Figure VI. 3.1 : Intersection of segments, which are not 
parallel to the axes 


where 0 is the angle between the X-axis and the normal to from 

the origin. Rotate the plane by -0, and denote the new quantities by 
primes, as shown in Figure VI. 3. 2. 



Figure VI. 3. 2 : Figure VI. 3. 1 after the rotation 
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The relation between P* and P is: 


(x 1 , y') = (x, y) 


COS 0 


sin 0 


■sin 0 
cos 0 


Pq is found easily as PJP^ is parallel to the Y-axis: 

P 0 = ^ x 0’ y 0^ where x o “ X 1 and y 0 = S{x A _x 0’ y A’ x B" x 0’ y B } 

Note that 6 multiplications are saved, (out of the 20) as y j, x^, y^ 
are not needed. 


From triangle similarity we get: 
'A 


x* - x* 

A x 0 


y A - Y'o 


l p o p il 


x 


B 


k 0 


y 


B 


'0 


P'P' 

1 0 B 


and 

X A " x 0 = y A : y 0 _ 1 P 0 P a1 

X B “ x 0 y B y 0 |P Q P B | 


where |P.P. | denotes the length of the segment P.P.. As 


i J 


i J 


I p o p aI = I p o p aI 


and 


P ! P ! = P P 
'OB* 1 OB 


we get: 

X A ~ x 0 = y A ~ y 0 = X A ~ x 0 
X B " x 0 y B " y 0 X B " X 0 

By using only 6 multiplications we can get xj^, x^, x^ and then finding 

p o by: 

x 0 = S{x A " x 0> X A’ X B " x 0’ X B } 

and 

y 0 = S{x A ‘ x 0’ y A> X B " x 0> y B } • 


This eliminates the need for back-rotation. Note that both SI*s above 
have the same values in their X-registers, and can be combined to one 
unit, having 6 registers, and 3 adders. 
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Next, we want to show how we can do without sin 0 and cos 0. 


Consider P^P^: 


Ax - x^ “ x^ = -X sin 0 


) where X = tjPjP^I 


Ay ~ y £ ~ V i = X cos 0 

From x' = x cos 0 + y sin©, one gets: 

Xx* = x- Ay - y* Ax . 

The SI is homogeneous in the x coordinates, so we have: 


where: 


II 

0 

X 

S ^ Xx A " 

Xx 0 ’ 

X A’ 

Xx B ' 

’ Xx <0, 

X B ) 

II 

0 

S ^ Xx A ” 

Xx 0 ’ 

*A’ 

Xx B ' 


y B } 

Xx A 

= x A - Ay 

■ y A' 

■ Ax 




Xx B 

= x B ' Ay 

' y B 

Ax 




Xx o 

= Xj- Ay 

" Yx* 

Ax 

* 




The SI rejects the case in which the intersection is outside P A P^ , 

_ A Jd 7 

but does not reject if the intersection is outside P^P^. If one 
cares for this condition as well, he should check for it when the re¬ 
sult is announced. 

If many segments have to be intersected with P^P^, the problem 
is solved using only 4 parallel multiplications per segment, as 

Xx o = x i A y " yi Ax 

depends only on P^ and P^. 

Note that the above technique does not treat P^Pg and P^P^ 
symmetrically, and is oriented to consider P^P^ as a segment but 
P^P^ as a line which causes the need for the extra checking for 
Pq e (Pp P^). If time is critical, and if we wish to consider both 
P^Pg and P^P^ as segments, then we use 2 similar units which compute 
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in parallel both (xj^ - x”, x^ - x^) and (xj - x^, x£ - x^) 1 . The trivial 
rejection can be done simultaneously for both cases. 

Another trivial rejection can be used by checking, prior (or in 
parallel) to any other computation the condition: 


( X A > Xj). AND. (x B > Xj). AND. (x A > x 2 ). AND. (x B > x 2 ) 

. OR. (x A < x^). AND. (x B < x^). AND. (x^ < x 2 ). AND. (xg < x 2 ) 

. OR. (y A > y^). AND. (y B > y j). AND. (y A > y 2 ). AND. (y B > y 2 ) 

• OR- (y A < yj). AND. (y B < y^. AND. (y A < y 2 ). AND. (y fi < y 2 ) . 

If this condition holds, then obviously there is no intersection. In the 
example of Figure VI. 3. 3 we illustrate all the trivial rejections. 



Figure VI. 3. 3 : The segments {P A , Pg} are trivially rejected 


1 The doub le prime quantities are in a system which was rotated such 
that is parallel to the Y-axis. 
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(VI. 4) On String-Line Intersecting 


In many real time applications it is essential to find in a very 
short time, all of the intersections of a given line with some sets of 
contiguous segments. We call such contiguous segments set "strings”. 

We describe below a device which can perform this task. This 
device is initialized by loading information about the line, then it is 
given the data about the end-points of the segments, and it finds all 
the intersections of these segments with the line specified earlier. 

The time consumed by each point is between one multiplication-time 
plus two addition times to one divide-time . With today*s technology 
one divide time is about the time required for fetching information from 
memory. 

The operation of this device is based on the method described in 
(VI. 3). The device described here is able to handle some operations 
in parallel. All the time-consuming operations, (namely fetching the 
data, multiplication and interpolation) are organized such that each 
of them starts as soon as possible, taking advantage of cases in which 
the other operations are fast. Because these operations can work com¬ 
pletely in parallel the time required per point is the duration of the 
longest operations. 

Let the strings be represented in memory by a sequence of triples 

{b., x i ? y^} where (x.,y.) are the coordinates of P., the i~th point, and 

b. is a bit which is set to 'one 1 if there is a segment between P. , and 
1 l-l 

P and it is ’zero 1 otherwise. Consider the example in Figure VI. 4. L: 
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Figure VI. 4.1 : A line and a string 

These segments are stored as: 

(°, Xj, Vi)(1, * 2 , y 2 )(l, x 3 , y 3 )(l, x 4 , y 4 )(0, x 5 , y 5 )(l, x & , y 6 )( 1, x ? , y ? ) . 
The line P^P^ to be intersected is defined by 

C = My 2 * y^ 
s = X(x 1 - x 2 ) 

6 = -(Cx 1 + Sy^ 

where X is a normalization factor. These quantities are calculated, 
once per line, and loaded to the C, S and 6 registers. 

After the initialization of the device the first triple is given to 
it, which starts to process the data. When it is ready to accept the 
next triple it calls its input channel. When an intersection is found 
it calls the output channel to record the intersection point. The or¬ 
ganization of the hardware is oriented toward a high degree of parallel 
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processing. The data about each point has to pass 2 multiplications 
in parallel and 2 additions in sequence, before it is given to the SI, 
together with the information about the previous point. If the b-bit 
is found to be 'zero*, the SI rejects the segment, otherwise it starts 
its operation which might reject the segment if there is no inter¬ 
section, or find the intersection point. In this case, the SI announces 
success and calls the output channel to record the intersection. 

The hardware organization is shown in Figure VI. 4. 2 and its 
operation is described by the graph in Figure VI. 4. 3. 
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Figure VI.4.2 : The line/segments intersector. 

































































V I . t \. S : The operat i n 
of the 1Lne/segmen ts 
intcrscctor. 
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